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SECTION  1 
INTRODUCTION 


1.1  OVERVIEW 

The  Earth  Satellite  Program  (ESP)  is  an  orbital  analysis  software  tool  for  Earth-orbiting 
satellites.  The  main  features  of  the  program,  in  a  broad  sense,  are  coverage,  tracking,  and 
mapping.  More  specifically,  die  program  has  the  capability  to 

(a)  generate  the  ground  track  of  a  satellite,  that  is,  the  projection  of  the  satellite  on  the 
surface  of  the  Earth  as  the  satellite  moves  in  its  orbit  over  time. 

(b)  determine  outage  zones  over  a  region  of  the  Earth  for  a  constellation  of  up  to  200 
satellites  and  a  user-specified  minimum  elevation  angle.  An  outage  zone  is  a  set  of 
ground  site  locations  that  are  unable  to  see  a  user-specified  minimum  number  of 
satellites  in  the  constellation  over  some  time  interval. 

(c)  project  single  or  multiple  antenna  spot  beams  of  a  user-specified  beam  width  in  any 
direction  onto  the  Earth’s  surface  and  track  the  location  of  these  beams  over  time. 
The  option  to  generate  contours  of  equal  loss  for  a  single  beam  is  available. 

(d)  calculate  the  elevation  angle  and  range  from  ground  site  locations  to  a  satellite,  each 
as  a  function  of  time.  Range  rate  is  also  included. 

For  each  of  these  features,  output  is  provided  in  tabular  form.  In  addition,  items  (a)  through 
(c)  may  be  plotted  on  all  or  a  portion  of  a  world  map.  Satellite  orbits  are  propagated  using  a 
fourth-order  Runge-Kutta  integration  scheme.  The  program  runs  on  a  Macintosh  computer 
and  is  written  in  THINK  Pascal  [4, 5].  The  software  on  diskette  is  also  attached. 


1.2  PURPOSE 

This  effort  evolved  from  some  computer  programs  with  similar  features  that  were 
developed  in  the  1980s  by  MITRE  D090.  Previous  programs  were  limited  to  circular,  two- 
body  orbits,  assumed  a  radially  symmetric,  point  mass  Earth,  and  were  prompt-driven.  In 
contrast,  ESP  accommodates  any  elliptical  orbit  and  includes  the  perturbations  due  to  a  non- 
spherical  Earth,  atmospheric  drag,  and  solar  radiation  pressure.  It  is  also  simple  to  use;  being 
Macintosh-based,  it  adheres  to  the  standard  Macintosh  look-and-feel  by  implementing  dialog 
boxes,  pull-down  menus,  menu  commands,  icons,  scroll  bars,  etc.  Furthermore,  unlike  its 
predecessors,  its  mapping  capability  is  an  integral  part  of  the  program. 

ESP  was  developed  for  the  Milstar  program  to  support  system  design  trade-off  analyses 
for  both  ground  and  space  segments.  In  particular,  it  may  be  used  to  determine  those  satellite 
constellations  or  ground  site  terminal  antenna  designs,  or  some  combination  thereof,  that 
provide  adequate  or  even  optimal  coverage  of  the  Earth.  Because  it  is  easy  to  learn  and  use  and 
typically  generates  results  in  a  matter  of  minutes,  the  program  is  a  valuable  tool  for  providing 


1 


quick  turn-around  responses  to  questions  that  may  arise  from  the  System  Program  Office 
(SPO).  However,  nothing  in  the  program  restricts  it  to  Milstar  applications.  ESP  is  quite 
general  in  nature  and  has  successfully  been  used  to  support  a  number  of  programs  at  the 
MITRE  Corporation. 

ESP  2.0  is  an  upgrade  from  version  1.0  which  was  released  in  June  1991.  Major 
enhancements  include  spot  beam  projections,  the  capability  to  enter  larger  constellations  and  to 
require  multiple  satellites  to  be  in  view  of  any  ground  site  in  the  outage  feature,  and  a  more 
sophisticated  mapping  capability  that  includes  map  zooming. 


1.3  REPORT  ORGANIZATION 

This  paper  documents  the  capabilities,  usage,  and  mathematical  models  of  ESP  version 
2.0,  including  recent  enhancements  and  all  features  in  version  1.0.  It  is  organized  into  five 
major  sections.  An  overview  of  the  content  of  each  of  the  sections  that  follows  is  presented 
here. 

Section  2  serves  as  a  user’s  guide.  It  begins  with  a  discussion  of  some  examples  of 
specific  kinds  of  problems  that  may  be  analyzed  using  the  program.  A  detailed  description  of 
how  to  operate  the  program  then  follows.  This  section  assumes  a  general  familiarity  with 
Macintosh-based  applications  that  rely  on  a  mouse. 

The  analytical  foundations  are  provided  in  section  3.  For  the  user  who  is  unfamiliar  with 
the  six  classical  orbital  elements,  which  are  required  inputs  for  all  aspects  of  the  program,  the 
initial  subsections  are  prerequisite  reading.  A  detailed  description  of  the  equations  of  motion, 
perturbations,  and  integration  algorithm  is  included.  This  section  also  explains  how  the 
satellite  position  vector  is  used  to  generate  the  ground  track  and  perform  the  elevation  angle, 
range,  and  range  rate  calculations,  and  how  the  spot  beam  projections,  horizon,  and  loss 
contours  are  derived. 

Section  4  contains  a  discussion  of  computer  programming  considerations.  Since 
computers  can  represent  numbers  with  only  finite  precision,  errors  in  floating  point 
calculations  are  inevitable.  This  section  details  how  ESP  handles  the  problems  of  floating 
point  calculations.  It  also  describes  how  the  program  was  tested  and  validated. 

Finally,  some  example  runs  illustrating  the  program’s  main  features  are  provided  in  the 
appendix.  This  section  includes  a  brief  description  of  the  various  types  of  orbits  analyzed.  All 
three  types  of  program  outputs  —  echoed  input  data  and  tabular  and  graphical  results  —  are 
shown.  However,  since  the  non-graphical  output  is  generally  quite  lengthy,  only  a  sample  of 
this  type  of  output  is  given  for  each  feature. 
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SECTION  2 
USER’S  GUIDE 


2.1  OVERVIEW 

This  section  describes  how  to  use  ESP  2.0.  The  program  presents  the  standard  Macintosh 
user  interface.  For  a  user  familiar  with  operating  Macintosh  programs  in  general,  exploring  the 
program  should  not  be  difficult  Specifically,  the  user  should  know  how  to  click  and  drag 
items,  invoke  pull-down  menu  commands,  size  windows,  and  use  other  features  common  to 
Macintosh  programs. 


2.2  PROGRAM  REQUIREMENTS 

ESP  2.0  requires  a  Macintosh  II  with  four  megabytes  of  random  access  memory  (RAM), 
System  Software  version  6.0.5  or  later,  and  a  high  density  floppy  disk  drive.  The  specific 
features  required  by  the  program  are  an  MC6888 1/882  mathematics  coprocessor,  Color 
QuickDraw1,  and  a  large  monitor  (at  least  12  inches  if  monochrome,  13  inches  if  color),  all  of 
which  are  standard  on  most  Macintosh  Os.  If  any  of  these  features  are  not  present,  ESP  will 
inform  the  user  of  the  problem  and  then  halt.  A  color  monitor  is  recommended  but  not 
required.  The  files  required  to  run  ESP  2.0  are  ESP  itself,  ESPmap,  and  ESPdensity.  These 
files  are  on  the  diskette  provided  with  this  document. 

Certain  Macintosh  initialization  files,  or  INITs,  that  operate  in  the  background  may 
interfere  with  ESP.  If  the  user  experiences  problems  running  ESP,  the  INTTs  should  be 
removed  until  ESP  runs  normally.  This  is  done  by  systematically  taking  them  out  of  the 
system  folder  and  restarting  the  Macintosh.  In  System  7.0,  INTTs  are  referred  to  as  system 
extensions.  They  can  be  found  in  the  Extensions  folder  within  the  System  folder. 


2.3  SAMPLE  USAGE 

Before  exploring  ESP’s  various  menus  and  commands,  some  examples  of  the  program’s 
utility  and  application  are  presented.  The  intent  of  this  section  is  to  provide  the  reader  with  a 
flavor  of  the  various  problems  that  may  be  analyzed. 

The  satellite  ground  track  indicates  the  location  of  the  satellite  with  respect  to  the  geography 
of  the  rotating  Earth.  Knowledge  of  this  location  is  important  when  addressing  various 
coverage  issues.  By  inspection,  one  can  obtain  a  general  sense  as  to  what  regions  of  the  Earth 
are  likely  to  have  coverage.  The  ground  track  may  be  used  as  a  tool  in  the  design  of  satellite 
constellations  that  provide  continuous  world-wide  coverage.  By  examining  the  individual 
ground  traces  of  several  different  orbits,  one  may  gain  insight  into  the  types  of  orbit 
combinations  that  optimize  coverage.  (This  is  especially  useful  in  conjunction  with  the  outage 


^olor  QuickDraw  is  the  part  of  the  Macintosh  system  software  that  supports  color  graphics.  It  is  not  the 
same  as  a  color  monitor.  Color  QuickDraw  is  not  present  on  some  older  Macintosh  models. 
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zone  feature.)  A  typical  run  propagates  a  satellite  trajectory  over  a  24-hour  period.  Longer 
runs  are  permitted  by  the  program;  however,  no  studies  have  been  made  as  to  how  the 
accuracy  of  the  model  degrades  over  time.  For  such  longer  runs  there  are  more  appropriate 
tools  available,  such  as  ASAP  listed  in  the  reference  section. 

The  elevation  angle  calculation  may  be  used  to  determine  how  low  a  ground  site  terminal 
antenna  must  be  pointed  in  order  to  maintain  coverage.  In  general,  the  elevation  angle  from  a 
given  ground  site  to  a  satellite  will  vary  as  the  satellite  moves  in  its  orbit  The  elevation 
calculation  will  indicate  the  elevation  angle  at  user-specified  times  over  the  course  of  the  run. 
The  ground  site  antennae  should  then  be  designed  so  that  the  ground  site  is  capable  of 
contacting  the  satellite  over  the  range  of  angles  calculated.  Alternatively,  the  elevation  angle 
calculation  may  be  used  to  assess  sensitivities  of  orbital  parameters  to  a  minimum  elevation 
angle  constraint.  For  example,  by  calculating  the  elevation  angles  from  ground  sites  to  a 
satellite  in  geosynchronous2  orbit  with  various  inclinations,  the  program  can  be  used  to 
determine  the  maximum  tolerable  inclination  drift  as  a  function  of  minimum  elevation  angle. 

The  range,  which  is  the  distance  from  a  given  ground  site  to  the  satellite,  is  useful  for 
determining  the  delay  time  between  the  transmitting  and  receiving  of  a  signal.  The  range  can 
be  divided  by  the  speed  of  light  to  obtain  the  time  required  for  the  signal  to  travel.  The  time 
derivative  of  range,  or  range  rate,  is  useful  in  determining  the  extent  of  the  Doppler  shift  on 
signals.  With  this  information  available,  the  ground  site  terminals  can  then  be  designed  to  take 
the  changing  frequencies  into  account. 

The  outage  zones  calculation  is  an  important  tool  for  determining  what  satellite 
configurations  and  antenna  designs  and  placements  will  result  in  optimal  coverage.  Several 
satellite  orbits  may  be  entered  at  once,  whereupon  ESP  will  determine  those  ground  sites  that 
are  not  covered,  and  an  approximation  of  the  outage  time  at  each  site.  The  orbit  of  each 
satellite  can  be  adjusted  individually  so  that  gaps  in  coverage  may  be  filled.  The  number  of 
satellites,  minimum  antenna  elevation  angle,  and  minimum  number  of  satellites  in  view  can 
each  be  adjusted  individually  in  order  to  assess  the  sensitivity  of  coverage  to  the  constellation’s 
configuration.  These  kinds  of  analyses  can  assist  in  determining  the  trade-off  between  more 
satellites  and  more  powerful  or  complex  antennae. 

The  spot  beam  feature,  in  addition  to  simply  determining  beam  coverage,  can  be  used  to 
assess  the  sensitivity  of  beam  coverage  to  satellite  placement  and/or  beam  width.  Since  a  given 
beam  width  corresponds  (approximately)  to  some  elevation  angle,  this  feature  can  also  be  used 
to  generate  elevation  angle  contours.  In  addition,  contours  of  equal  decibel  loss  show  how 
signal  strength  decreases  as  a  function  of  distance  from  the  beam  center. 


2.4  USING  THE  PROGRAM 

The  user  is  presumed  to  have  a  basic  knowledge  of  the  six  elements  of  an  orbit  and  of 
elementary  perturbed  two-body  theory.  An  adequate  background  can  be  obtained  by  reading 
section  3  of  this  document. 


2“Geosynchronous”  is  defined  in  appendix  A. 
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As  with  any  Macintosh  program,  ESP  is  started  by  double-clicking  its  icon  (see  figure  2-1) 
or  by  selecting  it  and  choosing  Open  from  the  Finder’s  File  menu.  In  order  for  the  program  to 
run  properly,  the  files  ESPmap  and  ESPdensity  must  be  in  the  same  folder  as  ESP.  If  the 
support  files  are  not  present,  ESP  will  still  run.  If  ESPmap  cannot  be  found  or  loaded,  ground 
tracks,  outage  zones,  and  spot  beams  will  still  be  plotted,  but  on  a  blank  rectangle.  The  only 
visible  geographical  frame  of  reference  will  be  the  longitude  and  latitude  labels  along  the  edges 
of  the  plot.  If  ESPdensity  cannot  be  found  or  loaded,  ESP  will  compute  the  orbits  as  if  there 
were  no  atmosphere.  After  ESP  loads  the  support  files,  a  dialog  box,  shown  in  figure  2-2  and 
described  in  section  2.4.2  under  New...,  is  displayed  asking  the  user  to  choose  a  calculation 
type. 


=| _ j  : - 

ss  ESP  Document  Icons  1 

i  —  gjs 

7  items 

1 01 .3  MB  in  disk 

52.9  MB  available 

A 


m 

ESP 


Ground 


ESPmap 


It 

ESPdensity 


rack  Document  Elevation  Angles  Document 


Outage  Zones  Document  Spot  Beam  Document 


o 


o 


a 


Figure  2-1.  ESP  Icons 


Which  type  of  calculation  mould 

you  like  for  the  neiu  document? 

®  Ground  track 

(i  °k  i 

O  Eleuation  angles 

O  Outage  zones 

Cancel  j 

Q  Spot  beam  projections 

Figure  2-2.  New  Document  Dialog  Box 
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The  user  may  also  invoke  ESP  from  any  previously  saved  ESP  documents  (saving  is 
described  in  section  2.4.2,  under  Save  and  Save  As. . .).  This  is  accomplished  by  selecting  one 
or  more  of  these  documents  and  then  either  double-clicking  on  one  of  them  or  choosing  Open 
from  the  Finder’s  File  menu.  The  method  is  the  same  as  when  opening  any  other  Macintosh 
document.  In  this  case  the  dialog  box  in  figure  2-2  will  not  be  displayed,  nor  will  the  user  be 
prompted  to  enter  new  inputs. 

Within  ESP,  five  menus  are  displayed  across  the  top  of  the  screen.  The  first  menu  is  the 
usual  apple  menu.  The  second  menu  is  the  File  menu,  which  contains  commands  for  opening, 
closing,  and  printing  documents.  The  third  menu,  the  Edit  menu,  is  not  used  in  the  current 
version,  but  is  fully  available  for  use  by  desk  accessories.  The  next  menu  is  the  Options 
menu,  which  is  used  to  invoke  any  of  the  program’s  four  main  features;  it  contains  all 
commands  for  calculations  and  drawing  maps.  Finally,  the  Windows  menu  enables  the  user  to 
reactivate  document  windows  that  have  been  temporarily  closed  or  hidden  by  other  windows. 

Occasionally  the  user  may  request  an  operation,  such  as  New  or  Open,  that  would  result  in 
more  windows  than  ESP  can  handle.  (The  limit  is  currently  100  windows  at  a  time.)  If  this 
happens,  a  message  will  appear  stating  that  too  many  windows  are  open.  In  this  situation,  the 
user  must  permanently  close  some  other  window,  as  described  in  section  2.4.3. 

2.4.1  The  Apple  Menu 

About  ESP,., 

This  command  displays  the  version,  authors,  and  copyright  notice. 

Desk  Accessories 

Any  desk  accessory  installed  on  the  system  may  be  invoked  as  usual.  When  a  desk 
accessory  is  opened,  the  Edit  menu,  normally  dormant  in  ESP,  is  made  available  for  copying, 
cutting,  and  pasting.  Some  of  the  commands  in  the  File  menu,  which  are  inappropriate  for  a 
desk  accessory,  are  disabled.  If  the  Macintosh  is  running  MultiFinder  or  System  7.0,  the  desk 
accessory  will  temporarily  hide  ESP’s  menus  and  substitute  its  own. 

2.4.2  The  File  Menu 

New.,., 

The  New. . .  command  creates  a  new  untitled  window.  A  dialog  box  is  presented  for 
selecting  the  desired  type  of  calculation;  ground  track,  outage  zones,  elevation  angles,  or  spot 
beam  projections  (see  figure  2-2).  The  elevation  angles  option  includes  range  and  range  rate 
calculations.  After  a  calculation  type  is  selected,  two  dialog  boxes  are  displayed  in  succession 
requesting  inputs  for  the  orbital  calculations.  These  dialog  boxes  are  described  in  section 
2.4.5. 
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Open... 


The  Open. . .  command  presents  the  standard  dialog  box  for  opening  documents.  Only 
documents  created  by  ESP  will  be  listed.  A  particular  document  is  opened  by  either  selecting 
its  name  in  the  list  and  clicking  the  “Open”  button  or  by  double-clicking  its  name.  If  a 
document  was  created  by  an  earlier  version  of  ESP,  the  program  will  notify  the  user  that  if  the 
document  is  saved  again,  it  will  be  saved  in  version  2.0  format  The  version  number  of  a 
document  can  be  determined  from  the  Finder  in  the  usual  manner  with  the  Get  Info  command 
under  the  File  menu. 

Close 

This  command  closes  the  active  window.  If  the  window  is  the  only  window  open  to 
particular  document  and  the  document  has  not  been  saved  since  its  last  calculation,  the  user  will 
be  given  an  opportunity  to  save  the  document. 

Save 

This  command  saves  the  active  window  and  its  inputs.  If  the  document  represented  by  the 
window  has  both  a  map  and  tabular  output,  both  representations  of  the  data  are  automatically 
saved  as  one  file  in  addition  to  the  input  dialog  boxes.  If  the  document  does  not  yet  have  a 
title,  this  command  is  equivalent  to  the  Save  As. . .  command. 

SaveAs... 

This  command  saves  the  active  window  under  a  new  title.  A  dialog  box  will  appear  for 
entering  the  desired  name.  The  new  title  will  be  reflected  in  the  document’s  window  or 
windows.  If  a  document  already  exists  on  the  disk  with  the  same  name,  the  user  will  be  given 
the  choice  of  either  replacing  the  existing  document  or  choosing  a  different  title. 

Revert  to  Saved,,, 

The  Revert  command  has  not  been  implemented  in  this  release  of  ESP.  One  can,  however, 
retrieve  the  saved  version  of  a  document  by  doing  the  following.  First,  all  windows  associated 
with  the  document  must  be  closed.  A  message  will  then  appear  asking  the  user  whether  or  not 
to  save  the  document.  The  user  should  click  the  “No”  button,  then  choose  Open  from  the  File 
menu,  and  select  the  document. 

Page  Setup. . . 

This  command  presents  a  dialog  box  that  allows  the  user  to  change  the  printer  settings. 

The  dialog  box  will  be  different  for  different  printers.  Example  settings  that  can  be  changed 
are  reduction/enlargement,  page  orientation,  and  printer  effects. 

faint.-,. 

This  command  enables  the  user  to  print  the  active  document’s  calculation  results,  and  the 
map,  if  the  document  contains  ground  track,  outage  zone,  or  spot  beam  calculations.  If  a  map 
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is  available,  a  dialog  box  will  be  presented  asking  which  portions  of  the  output  are  to  be 
printed  (see  figure  2-3).  The  output  is  divided  into  three  categories:  the  echoed  input,  the 
tabular  output,  and  the  map.  This  feature  was  designed  so  that  the  user  may  print  the  inputs 
and  the  textual  output  on  a  laser  or  impact  printer  and  the  map  on  a  color  printer.  It  is  also 
useful  if  the  user  desires  only  part  of  the  output.  If  the  computer  is  on  a  network  with  more 
than  one  printer,  the  Page  Setup  command  should  be  executed  after  switching  to  a  different 
printer.  Once  the  output  types  to  be  printed  are  selected  (or  immediately  after  selecting  Print, 
in  the  case  of  an  elevation  angles  document),  the  standard  print  dialog  box  is  presented  and 
asks  for  page  numbers  and  other  information  specific  to  the  job.  Print  spooling  is  supported. 
For  legibility,  the  printout  uses  a  different  character  font  from  that  on  the  screen. 


Which  portions  of  “Untitled  1”  do  you 
want  printed? 


^  Echo  the  inputs 
^  Print  tabular  output 
£3  Draw  the  map 

Cancel  j 


c 


Print 


3 


Figure  2-3.  Print  Options  Dialog  Box 


Quit  ESP 


This  command  terminates  execution  of  ESP.  If  any  documents  have  not  been  saved,  the 
user  is  given  the  opportunity  to  save  them. 

2.4.3  Window  Manipulation 

This  subsection  describes  properties  of  ESP  windows  that  are  not  present  in  all  Macintosh 
programs. 

Although  the  height  of  a  text  window  may  be  changed  by  dragging  the  size  icon  in  the 
lower  right  comer  of  the  window,  the  width  cannot  be  changed.  (Both  the  width  and  the 
height  of  a  map  window  can  be  adjusted,  as  in  most  Macintosh  applications.) 


8 


Outage  zone  and  spot  beam  documents  have  a  horizontal  scroll  bar  at  the  bottom  of  their 
text  windows.  This  scroll  bar  does  not  scroll  between  columns  but  is  used  instead  to  switch 
between  outage  zones  or  individual  beam  projections. 

Holding  down  the  command  key  (8€)  or  the  option  key  while  clicking  the  close  box  in  the 
upper  left  comer  of  a  window  hides  the  window  temporarily.  The  window  is  still  in  memory, 
however,  and  may  be  recalled  using  the  Windows  menu.  Clicking  the  close  box  without 
holding  down  either  of  these  keys  closes  the  window  permanently.  Note  that  this  behavior  is 
different  from  that  of  ESP  1.0,  where  clicking  the  close  box  without  the  command  or  option 
key  hid  the  window  temporarily. 

If  a  message  appears  asking  the  user  to  close  some  windows,  the  user  must  close  a 
window  permanently,  i.e.,  without  using  the  command  or  option  keys. 

2.4.4  The  Edit  Menu 

At  this  time,  ESP  does  not  use  the  editing  commands  (Undo,  Cut,  Copy,  Paste,  and 
Clear).  However,  these  commands  arc  made  available  when  a  desk  accessory  is  opened. 

2.4.5  The  Options  Menu 

The  first  four  commands  in  this  menu  are  Calculate  Ground  Track. . .,  Calculate  Elevation 
Angles. . .,  Calculate  Outage  Zones. . .,  and  Calculate  Spot  Beam  Projections....  After  choosing 
one  of  these  commands,  the  user  goes  through  the  same  process  for  each.  Two  dialog  boxes 
appear  in  succession,  in  which  the  user  fills  in  information  about  the  satellite  or  constellation 
and  the  particular  run.  Not  all  the  input  fields  are  needed  for  the  ground  track,  elevation  angle, 
or  spot  beam  calculations,  so  certain  fields  are  grayed  out  in  these  cases.  Any  information 
entered  into  these  fields  will  be  ignored.  The  dialog  boxes  are  shown  in  figures  2-4  through 
2-11  with  their  default  entries. 

The  first  dialog  box  solicits  information  about  the  satellite  or  constellation.  The  scroll  bar 
and  the  field  for  the  number  of  satellites  are  used  in  defining  constellations  of  satellites.  Since 
constellations  are  used  only  for  the  outage  zone  calculations,  the  discussion  of  these  fields  will 
be  deferred  to  the  section  on  outage  zones.  The  next  seven  fields  are  for  inputting  the  six 
classical  orbital  elements  of  the  satellite.  These  elements  are  explained  in  section  3.2.1.  As 
indicated  in  the  dialog  box,  all  angles  are  in  degrees,  and  distance  is  in  kilometers.  The 
inclination  must  be  between  0*  and  180*  inclusive.  The  longitude  of  the  node  and  the  argument 
of  perigee  must  be  greater  than  or  equal  to  0*  and  less  than  360*.  Dates  in  the  dialog  boxes 
must  be  in  the  form  YYYYMMDD  (e.g.  19930704  for  4  July  1993)  and  times  must  be  in  the 
form  HHMMSS  using  twenty-four  hour  time  (e.g.,  130000  for  1:00  PM). 

Note  that  the  longitude  of  the  node,  referred  to  above,  is  measured  relative  to  the  vernal 
equinox  and  is  not  the  same  as  the  more  familiar  longitude  system,  which  is  based  on  the 
Greenwich,  England  prime  meridian.  The  vernal  equinox  is  based  on  a  coordinate  system  that 
is  fixed  with  respect  to  the  stars,  whereas  the  prime  meridian  rotates  with  the  Earth.  Hence, 
the  subsatellite  longitude,  which  is  relative  to  the  prime  meridian,  is  in  general  often  not  equal 
to  the  longitude  of  the  node.  For  example,  for  the  default  date  and  time  inputs  the  difference 
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?%u  Ground  Track 

ir  Please  enter  the  satellite  Information: 


Satellite*  I : 


'[tyl  Number  of  satellites: 


Semimajor  am*  (a,  km): 
Eccentricity  (eh 
Inclination  (l,  deg.): 

Longitude  of  the  node  (0,  deg.): 
argument  of  perigee  (m,  dag.): 
Date  and  time  of  perigee  (t): 


0.0 


0.0 


0.0 


0.0 


19910101 


c 


Calculate 


Cancel 


Enter  more 


Input...  ) 


000000 


19910101 


Oete  and  time  of  osculating  values: 

Mass  of  satellite  (kg): 

Effectlve  area  for  drag  (sq.  m): 

Coefficient  of  drag: 

Effective  area  for  solar  radiation  pressure  (sq.  m): 
Rbsorptlon/reflectlulty  constant: 


000000 


2.0 


toto/tt nw  cofttbfruttonf  should 
h*v*  ttve  format  YYYVMMDO 
rtHISS 


1.5 


Figure  2-4.  First  Ground  Track  Dialog  Box 


Ground  Track 

££\  Please  enter  the  following  additional  information: 


Latitude  of  firsl  site  (deg.i: 
lalllude  of  last  site  ideg.): 
Latitude  slop  size  (deg.i: 

Longitude  of  firsl  site  (deg.>: 
longilude  of  last  she  ideg.): 
Longitude  slep  size  (deg.i: 

Starting  date  and  time: 
Ending  date  and  time: 

Time  step,  In  hours: 

Minimum  antenna  angle: 
Hours  in  oath  outage  zone: 
Minimum  no.  of  satellites: 


Figure  2-5.  Second  Ground  Track  Dialog  Box 
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Elevation  Angles 

Please  enter  the  satellite  information: 


Satellite  *  I:  K>|  |iy|  Number  o 

Semimajor  ohIs  (a,  km):  G3BXH 

Eccentricity  (e):  (LO _  ~ 

Inclination  (I,  deg.):  0.0  Cl 

Longitude  of  the  node  (ft,  deg.):  OJ) _ 

Argument  of  perigee  (a,  deg.):  IM) _ _ 

Date  and  time  of  perigee  (t):  19910101  000000 


Number  of  satellites:  i 

1 _ 

1  [  Calculate  ~ 

[  Cancel 

_  fToteTmorTlnpurT 


Date  and  time  of  osculating  values:  19910101  000000 


Mass  of  satellite  (kg):  _  , 

Effective  area  for  drag  (sq.  m):  1 

Coefficient  of  drag:  2^0 _ 

Effective  area  for  solar  radiation  pressure  (sq.  m): 
Absorption/reflectivity  constant: 


Date/time  combinations  should 
have  the  format  YYYYMtDO 
NHhWSS. 


Figure  2-6.  First  Elevation  Angles  Dialog  Box 


Elevation  Angles 

Please  enter  the  following  additional  information: 


Latitude  of  first  site  (deg.):  |QjH 
Latitude  of  last  site  (deg.):  50.0 

latitude  step  size  (deg.):  5.0 

Longitude  of  first  site  (deg.):  -1 25.0 
Longitude  of  last  site  (deg.):  -65.0 
Longitude  step  size  (deg.):  5.0 

Starting  date  and  time:  199I0H 

Ending  date  and  time:  19910L 

Time  step,  In  hours:  0.5 

Minimum  antenna  angle:  [ 

Hours  in  eac  h  outage  zone: 

Minimum  no.  of  satellites: 


Calculate 


Previous  window., 


19910101 

000000 

19910102 

000000 

Negative  values  for  longitude 
and  latitude  Indicate  west  or 
south.  Date /lime  combinations 
should  have  the  format 
YYVYMMDD  HHMMSS. 


Figure  2-7.  Second  Elevation  Angles  Dialog  Box 
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V  Outage  Zones 

Li  Please  enter  the  following  additional  information: 


Latitude  of  first  site  (deg.): 

Latitude  of  last  site  (deg.):  50.0 

latitude  step  size  (deg.):  5.0 

Longitude  of  first  site  (deg.):  -125.0 
Longitude  of  last  site  (deg.):  -65.0 
Longitude  step  size  (deg.):  5.0 

Starting  date  and  time:  199101 

Ending  date  and  time:  199101 

Time  step,  in  hours:  0.5 

Minimum  antenna  angle:  10.0 

Hours  in  each  outage  zone:  3.0 
Minimum  no.  of  satellites:  I 


Calculate 

_ Cancel 

Previous  window.. 


19910101 

000000 

19910102 

000000 

Negative  vjlwi  for  longitude 
end  latitude  indicate  west  or 
south.  Dote /time  combinations 
should  hove  the  format 
YYYYMMDO  HHMMSS 


Figure  2-9.  Second  Outage  Zones  Dialog  Box 


Spat  Baam  Projections 
Please  anter  tba  satellite  Inrarmatlaa: 


Satellite  *1:  Ki»|  iCh  Number  o 

Semimajor  auto  (a,  km):  DDDH  r~ 

Eccaatrtclty  (ah  0.0 

Incllnatlen  (I,  Bag.):  li _  C 

Longitude  of  the  aeda  (O,  Bag.):  0^0 _  r 

argument  af  perigee  im,  dag.h  0J _ _ 

■ata  aad  ttma  af  partgaa  (r):  19910101  000000 


Number  af  to  tr  Hites:  |  | 


Calculata 

Cancal 

Enter  mora  Input.. 


Data  aad  tima  af  asculatlng  values:  [19910101  |000000| 


Mata  af  aatalllta  (kgh  _  , 

iffactiwa  araa  far  drag  (ag.  m):  _  1 

Coefficient  af  drag:  2^0 _ 

Effactiva  araa  far  aalar  radiatlan  praaaura  (aq.  m): 
Obsorptloa/roflectUilty  caaatant: 


Otto  Atom  MiMuHm  AnM 
key*  «w  (nrt  YYYYPPCO 


Figure  2-10.  First  Spot  Beam  Projections  Dialog  Box 


^  Pla 


Spat  Baam  Prajactiona 
’  the  fallowing  additional  Information: 


Baam  #1:  £J _ l$j  Numbar  of  beama:  [|1| 

Baam  longltuda  (dag.):  -99.41 

Baam  latitude  (deg.h  *0 _  , - cancel - ) 

Beam  width  (dag.):  4.0  _ _ _ 

_  [  Previous  window...  ] 

Starting  data  and  time:  19910101  000000 
Ending  data  and  tima:  19910101 loooooo] 

Time  step.  In  hours:  0.5  j^j^y/^J^y****"** 


Baam  longltuda  (dog.h  -99.41 
Beam  latitude  (deg.h  0.0 
Boom  width  (dag.):  4.0 


Time  atop,  In  hours:  |0.5 

□  Show  contours  of  oqual  loss 
Contour  step  sl?a:  [  ~  1 


O  Dogroot 
ODa«  Ibels 


•f  total  toes  «rt  tvaUto  tnly 
If  Start  U  tontto  m  toatnt. 


Figure  2-11.  Second  Spot  Beam  Projections  Dialog  Box 


between  the  two  values  is  99\41 .  A  program  that  calculates  the  longitude  of  the  ascending 
node  for  a  geosynchronous  satellite  centered  on  a  particular  subsatellite  longitude  is  supplied 
on  the  diskette  and  is  described  in  section  2.6. 

ESP  uses  ephemeris  time  (ET)  rather  than  universal  time  (UT).  Ephemeris  time  is  a 
uniform  measure  defined  by  the  orbital  motions  of  the  planets.  Universal  time  (sometimes 
referred  to  as  Zulu  or  Greenwich  mean  time),  on  the  other  hand,  is  based  on  the  diurnal 
rotation  of  the  Earth  and  the  motion  of  the  Earth  in  its  orbit  about  the  sun.  Differences  between 
the  two  measures  at  a  given  instant  are  slight,  generally  no  more  than  a  few  minutes  in  the 
twentieth  century.  Reference  16  provides  the  following  approximation  for  the  difference  (in 
minutes)  for  dates  beyond  1987: 

AT  =  ET  -  UT  =  0.41  + 1.20537'  +  0. 4992T2 

where  T  represents  the  time  elapsed  in  centuries  from  noon  on  31  December  1899  (Julian  date 
2415020.0)  to  present.  Exact  values  for  past  years  and  a  projected  value  for  1993  can  be 
obtained  from  The  Astronomical  Almanac  for  the  Year  1993  [17].  This  almanac  is  updated 
annually. 

The  default  semi-major  axis  is  42,163  km,  which  is  the  radius  of  a  geosynchronous  orbit. 
(Geosynchronous  orbits  are  defined  in  section  A.l.)  The  eccentricity,  inclination,  longitude  of 
the  node,  and  argument  of  perigee  are  all  initially  set  to  zero.  The  default  date  and  time  of 
perigee  is  January  1, 1991  at  midnight. 

The  date  and  time  of  the  osculating  values,  or  epoch  time,  is  the  time  at  which  the  six 
orbital  elements  are  known.  It  is,  in  effect,  the  time  of  the  initial  conditions.  Due  to 
perturbations  on  the  orbit  of  the  satellite,  the  satellite’s  six  orbital  elements  are  functions  of 
time,  so  this  information  is  important.  The  date  and  time  of  the  osculating  values  must  be 
before  or  at  the  date  and  time  of  perigee.  The  default  is  January  1,  1991  at  midnight. 

The  remaining  fields  are  for  inputting  information  that  is  required  for  atmospheric  drag  and 
solar  radiation  pressure  calculations.  To  include  these  effects,  positive  values  for  the  mass  of 
the  satellite  (in  kilograms)  and  the  effective  areas  (in  square  meters)  must  be  entered.  To  make 
a  run  without  one  or  both  of  these  effects,  the  appropriate  field(s)  should  be  filled  with  a  zero 
or  left  blank.  Both  of  these  effects  may  also  be  ignored  by  leaving  the  mass  field  blank.  If  an 
area  field  is  omitted  or  set  to  zero,  “N/A”  will  appear  on  the  printout  for  the  respective  field. 
“N/A”  will  also  be  printed  for  both  areas  and  the  mass  if  the  mass  is  omitted.  The  coefficient 
of  drag  has  a  default  value  of  2.0.  The  absorption/reflectivity  constant  ranges  from  0.0  to  2.0. 
It  assumes  a  value  less  than  1.0  for  varying  degrees  of  transparency,  1.0  for  a  total Jy  absorbent 
surface  (i.e.,  a  black  surface),  and  2.0  for  a  totally  reflective  surface.  An  example  is 
aluminum,  which  has  a  value  of  1.96.  The  program  uses  1.5  as  a  default. 

After  entering  the  satellite  information,  the  user  should  press  the  button  marked  “Enter 
more  input. . .”  ’Hie  second  dialog  box  will  appear.  This  dialog  box  concerns  itself  with 
aspects  of  the  particular  run.  Since  it  is  similar  for  the  ground  track,  outage  zone,  and 
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elevation  angle  calculations  it  is  described  here.  The  second  dialog  box  for  the  spot  beam 
calculation  is  very  different;  it  is  described  in  the  section  on  spot  beam  projections. 

For  the  elevation  and  outage  zones  calculations  the  first  six  fields  of  this  second  dialog  box 
determine  a  quadrille  grid  of  ground  sites  to  be  included  in  the  calculation.  The  longitudes  of 
the  first  and  last  sites  must  be  between  -180’  and  180*  inclusive,  where  the  positive  direction  is 
east  of  Greenwich  and  the  negative  direction  is  west  of  Greenwich.  The  latitude  entries  must 
be  between  -90*  and  90*  inclusive,  where  the  positive  direction  is  north  of  the  equator  and  the 
negative  direction  is  south  of  the  equator.  The  step  size  in  both  cases  must  be  positive.  The 
ending  longitude  and  latitude  are  included  in  the  grid  even  if  the  entire  range  in  either  case  is 
not  an  integral  multiple  of  the  step  size.  The  default  entries  form  a  grid  over  the  forty-eight 
contiguous  United  States  (CONUS)  with  sites  at  five-degree  intervals.  A  run  for  a  single  site 
can  be  generated  by  setting  the  ending  latitude  and  longitude  equal  to  the  starting  latitude  and 
longitude,  respectively. 

ESP  uses  geocentric  latitude,  which  is  the  angle  between  the  equatorial  plane  and  the  line 
passing  through  the  ground  site  and  the  center  of  the  Earth.  Geodetic  latitude,  which  is  often  the 
basis  for  maps  [2],  is  the  angle  between  the  equatorial  plane  and  the  line  perpendicular  to  the 
surface  of  the  Earth  at  the  ground  site  (see  figure  2-12).  Due  to  the  Earth’s  oblateness,  the  two 
latitude  values  of  a  given  ground  site  are  generally  different.  The  maximum  difference  between 
the  two  values  is  approximately  0.19  degrees,  which  occurs  at  -45*  and  45*  geodetic  latitude  and 
corresponds  to  a  difference  of  about  13  miles.  As  the  latitude  approaches  the  equator  or  the 
poles,  this  difference  approaches  zero.  Geodetic  latitude  (L)  can  be  easily  converted  to 
geocentric  latitude  (L‘)  using  the  following  formula: 

L'  =  arctan  (0.9933056 tan L) 


Figure  2-12.  Geodetic  (L)  Versus  Geocentric  (Z/)  Latitude 

The  next  three  fields  are  the  starting  date  and  time,  the  ending  date  and  time,  and  the  time 
step.  The  starting  and  ending  dates  and  times  are  given  in  the  format  described  above.  The 
time  step  is  given  in  hours.  The  ending  time  is  included  in  the  calculation  even  if  the  entire 
time  of  the  run  is  not  a  multiple  of  the  time  step.  Example  time  step  values  are  given  in  the 
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appendix.  It  should  be  noted  that  the  time  step  field  is  used  for  the  output,  not  for  the  orbital 
calculations.  The  calculations  have  their  own  internal  time  step  of  approximately  one  minute. 

The  last  three  fields,  which  are  used  only  for  outage  zone  calculations,  are  described  in  the 
section  on  outage  zones. 

The  user  may  return  to  the  first  dialog  box,  if  desired,  by  pressing  the  button  marked 
“Previous  window. . .”  The  user  may  switch  back  and  forth  between  the  two  dialog  boxes  as 
often  as  needed.  Once  satisfied  with  the  inputs,  the  user  should  press  the  “Calculate”  button. 
After  a  few  moments  during  which  the  satellite’s  orbit  is  calculated  (or  orbits,  if  there  are  more 
than  one  satellite),  the  results  will  be  displayed  in  the  text  window.  The  outputs  for  each  type 
of  calculation  are  described  in  their  respective  sections.  If  a  map  window  is  already  displayed, 
it  will  be  updated. 

Calculate  Ground  Track. . . 

This  option  calculates  an  individual  satellite’s  ground  track  and  altitude  as  a  function  of 
time.  It  is  only  available  when  the  active  window  is  a  ground  track  window.  The  ground  track 
dialog  boxes  are  shown  in  figures  2-4  and  2-5.  Since  only  one  satellite  is  accommodated,  the 
scroll  bar  and  the  field  for  the  number  of  satellites  are  grayed  out  in  the  first  dialog  box.  In  the 
second  dialog  box,  only  the  starting  date  and  time,  the  ending  date  and  time,  and  the  time  step 
are  required. 

The  ground  track  output  consists  of  four  columns.  These  columns  contain  the  date  and 
time,  the  subsatellite  longitude,  the  subsatellite  latitude,  and  the  satellite’s  altitude,  respectively. 
If  at  some  time  the  satellite  crashes  into  the  Earth,  the  user  will  be  alerted  with  a  dialog  box. 
After  the  user  dismisses  the  dialog  box,  the  output  line  at  that  point  will  be  displayed  in  red  (if 
a  color  monitor  is  used)  and  the  calculation  will  halt. 

Calculate  Elevation  Angles, , , 

For  a  given  set  of  ground  site  locations,  this  option  calculates  the  elevation  angles  and 
ranges  to  a  satellite  over  time.  Range  rate  as  a  function  of  time  is  also  included.  This  option  is 
available  if  the  active  window  is  an  elevation  window.  The  elevation  dialog  boxes  are  shown 
in  figures  2-6  and  2-7.  As  is  the  case  with  the  ground  track  dialog,  the  scroll  bar  and  the  field 
for  the  number  of  satellites  are  grayed  out  in  the  first  dialog  box.  The  information  to  be  entered 
is  the  same  as  for  the  ground  track.  In  the  second  dialog  box,  all  fields  are  used  except  the  last 
three. 

The  elevation  output  consists  of  six  columns.  Each  line  of  the  output  will  show  a  date  and 
time,  the  site  longitude,  the  site  latitude,  the  elevation  angle,  the  range,  and  the  range  rate.  If 
an  elevation  angle  is  negative  and  a  color  monitor  is  used,  the  angle  will  be  displayed  in  blue. 

.Calcuiats.Outagc  Zones, , , 

This  option  determines  the  coverage  provided  by  a  constellation  of  satellites  over  time.  It  is 
available  when  the  current  window  is  an  outage  zone  window.  The  coverage  may  apply  to  the 
entire  Earth’s  surface  or  to  a  portion  that  can  be  defined  by  a  grid  of  latitude  and  longitude 
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coordinates.  The  first  dialog  box  is  shown  in  figure  2-8.  In  addition  to  the  fields  described 
above,  the  first  two  fields  are  made  available.  The  number  of  satellites  field  is  for  setting  the 
size  of  the  constellation.  Up  to  200  satellites  may  be  entered.  Since  each  satellite  has  its  own 
orbital  elements  and  characteristics,  this  information  must  be  entered  separately  for  each. 
However,  the  information  entered  for  the  first  satellite  is  inherited  by  all  subsequent  satellites  in 
the  constellation.  The  scroll  bar  allows  the  user  to  switch  to  the  remaining  satellites  and  edit 
their  characteristics.  If  the  number  of  satellites  is  increased  in  a  document  that  was  read  from 
the  disk,  the  new  satellites  will  inherit  the  characteristics  of  the  last  satellite  in  the  original 
constellation. 

The  second  dialog  box,  shown  in  figure  2-9,  uses  all  the  input  fields.  The  minimum 
antenna  angle  determines  how  low  a  satellite  can  be  with  respect  to  the  horizon  while  still 
maintaining  coverage.  When  the  elevation  angle  from  a  given  ground  site  to  a  satellite  is  below 
this  minimum,  the  satellite  cannot  communicate  with  the  ground  site  and  will  be  considered 
blacked  out.  Only  a  single  value  can  be  entered;  this  will  be  used  for  all  satellites  and  ground 
site  locations. 

An  outage  zone  is  a  set  of  ground  sites  that  lack  coverage  from  a  minimum  number  of 
satellites  in  the  constellation  for  certain  amounts  of  time.  The  hours  in  each  outage  zone  field 
determines  the  size  of  the  interval  for  each  time  length  category.  For  example,  if  there  were 
three  hours  in  each  outage  zone,  then  zone  one  would  contain  all  sites  that  are  out  for  up  to 
three  hours,  zone  two  would  contain  all  sites  that  are  out  between  three  and  six  hours,  etc.  The 
program  can  accommodate  a  maximum  of  six  zones;  therefore  the  last  zone  often  contains  sites 
that  are  out  over  a  longer  time  interval.  Assuming  the  above  example  is  over  a  twenty-four 
hour  time  period,  zone  six  would  contain  those  sites  that  are  out  between  fifteen  and  twenty- 
four  hours.  Each  range  is  inclusive  of  the  upper  bound,  but  exclusive  of  the  lower  bound. 

For  example,  a  site  with  no  coverage  for  exactly  three  hours  would  be  in  zone  one,  not  zone 
two.  A  site  that  is  covered  for  the  entire  simulation  run  is  not  considered  to  be  in  any  zone. 

The  coverage  at  each  site  is  checked  at  the  user-defined  time  intervals.  In  other  words,  if 
the  user  specified  a  step  size  of  a  half  hour,  the  coverage  would  be  checked  every  half  hour. 
Thus,  the  outage  time  calculated  should  be  considered  an  approximation.  If  a  site  is  covered  at 
both  ends  of  a  time  interval,  it  is  assumed  to  be  covered  for  die  entire  interval.  If  a  site  cannot 
contact  at  least  the  minimum  number  of  satellites  at  both  ends  of  a  time  interval,  it  is  assumed 
to  bt  clacked  out  for  the  entire  interval.  If  a  site  is  covered  at  only  one  end  of  a  time  interval,  it 
is  assumed  to  be  covered  for  half  the  interval.  The  total  outage  time  for  a  site  is  thus  an  integral 
multiple  of  one-half  the  time  step. 

The  last  input  is  the  minimum  number  of  satellites  required  to  be  in  view  (above  the 
minimum  elevation  angle)  in  order  for  a  site  to  be  considered  covered. 

When  the  calculation  is  completed,  the  text  window  will  contain  three  columns  listing  the 
site  longitude  and  latitude  and  an  approximate  total  outage  time  for  each  sit;  th?  was  at  least 
partially  out.  The  output  itself  is  organized  by  outage  zone.  The  user  can  switch  between 
zones  by  manipulating  the  horizontal  scroll  bar  at  the  bottom  of  the  window. 
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This  option  generates  projections  of  satellite  spot  beams  onto  the  Earth  and  optionally 
generates  contours  of  equal  decibel  (dB)  loss.  The  first  dialog  box  is  the  same  as  for  the 
ground  track  and  elevation  angle  documents.  Only  one  satellite  may  be  entered  per  run. 

The  second  dialog  box  is  completely  different  from  that  of  the  other  calculation  types  (see 
figure  2-11).  The  user  may  specify  more  than  one  beam  from  the  satellite,  using  the  input  field 
labeled  “Number  of  Beams.”  The  scroll  bar  is  used  to  switch  between  beam  inputs.  The  beam 
longitude  and  latitude  specify  the  point  on  the  Earth  at  which  the  center  of  each  beam  is  aimed. 
The  beam  width  is  the  angle  that  subtends  the  beam  cone.  By  default,  there  is  one  beam  aimed 
approximately  directly  below  the  satellite,  with  a  width  of  four  degrees. 

The  starting  date  and  time,  the  ending  date  and  time,  and  the  time  step  specify  the  times  at 
which  spot  beam  calculations  are  made.  These  fields  are  the  same  as  those  so  labeled  in  the 
dialog  boxes  for  the  other  document  types. 

When  the  box  labeled  “Show  contours  of  equal  loss”  is  checked,  the  program  will 
generate  loss  contours  corresponding  to  a  single  projection  at  a  single  instant  in  time.  This 
feature  is  available  when  the  beam  width  is  30’  or  less.  A  circular  aperture  uniform 
distribution  beam  pattern  function  is  used.  The  mathematical  description  is  provided  in  section 
3.10.  When  this  option  is  selected,  certain  fields  are  grayed  out  (see  figure  2-13).  These  are 
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Figure  2-13.  Second  Spot  Beam  Projections  Dialog  Box  with  Contour  Input  Fields 


the  number  of  beams  field,  the  scroll  bar,  the  ending  date  and  time  fields,  and  the  time  step 
field.  The  contour  step  size  field  is  used  to  specify  the  intervals  for  generating  the  loss 
contours.  The  buttons  marked  “Degrees”  and  “Decibels”  define  its  unit.  By  default,  the 
“Degrees”  button  is  selected,  and  the  contour  step  size  is  set  to  one.  The  specified  beam  width 
is  assumed  to  correspond  to  -3  dB.  Additional  higher  loss  contours  are  generated  according  to 
the  step  size,  up  to  a  loss  value  of  -17.6  dB,  as  explained  in  section  3.10. 

To  show  contours  for  one  beam  from  a  document  that  contains  calculations  for  several 
beams,  the  user  should  scroll  to  the  desired  beam  before  checking  on  the  contour  loss  box. 

The  other  beams  will  then  be  discarded.  For  this  reason  it  is  recommended  that  the  user 
compute  the  contours  on  a  copy  of  the  original  multibeam  document.  (A  copy  can  be  easily 
made  with  the  Save  As. . .  command  from  the  File  menu  or  the  Duplicate  command  from  the 
Finder.) 

If  any  beam’s  pointing  location  is  not  visible  from  the  satellite  (i.e.,  falls  below  the 
horizon),  the  user  will  be  notified.  If  none  of  the  pointing  locations  are  visible  from  the 
satellite,  the  dialog  boxes  will  be  presented  to  the  user  for  editing. 

The  spot  beam  feature  displays  the  map  by  default  instead  of  the  text  window.  The  spot 
beam  projections  and/or  contours  are  displayed  as  solid  curves  on  the  map.  If  the  calculation 
was  done  over  time  and  the  user  has  a  color  monitor,  the  colors  of  the  curves  are  varied  to 
reflect  projections  at  different  times.  All  projections  at  a  given  time  are  displayed  in  the  same 
color,  though  the  color  associated  with  each  time  is  not  necessarily  unique.  When  plotting 
contours  of  equal  loss  each  contour  is  labeled  with  the  decibel  loss. 

If  any  beam  is  projected  partially  beyond  the  horizon  and  the  calculation  is  done  for  only 
one  instant  in  time,  the  satellite’s  horizon  contour  is  also  displayed.  The  horizon  is  represented 
as  a  dashed  blue  line. 

The  text  window  is  available  as  usual  via  the  Show  Text  Output  command  in  the  Options 
menu  (described  later  in  this  section).  The  text  window  displays  points  along  each  spot  beam 
projection,  and  points  along  the  horizon  or  loss  contours  if  applicable.  Other  information, 
such  as  beam  width  and  pointing  location,  is  displayed  above  the  column  headers. 

Draw  a  Map 

This  command  plots  the  data  from  a  text  window  onto  a  map  of  the  Earth.  The  command 
is  available  if  the  active  window  is  a  ground  track,  outage  zone,  or  spot  beam  text  window  and 
a  calculation  has  already  been  performed  for  that  window.  No  map  is  drawn  for  elevation 
angle  calculations.  If  the  map  has  already  been  plotted,  this  command  simply  brings  the  map 
window  to  the  front  of  the  screen. 

The  map  data  consists  of  approximately  200,000  points  taken  from  the  World  Data  Bank  II 
from  National  Technical  Information  Service.  The  data  includes  U.  S.  state  and  international 
boundaries,  current  as  of  the  breakup  of  the  Soviet  Union  (January  1992).  Germany,  Yemen, 
the  former  Czechoslovakia,  and  the  former  Yugoslavia  are  each  represented  as  a  single 
country.  The  map  data  has  a  resolution  of  one  minute  of  arc. 
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The  map  projection  used  is  equidistant  cylindrical,  a  projection  commonly  used  in  scientific 
plotting.  In  this  projection  both  the  meridians  of  longitude  and  the  parallels  of  latitude  are 
equally  spaced.  The  poles  are  represented  as  the  horizontal  lines  at  the  top  and  bottom  edges  of 
the  map.  The  longitude  and  latitude  scales  are  labeled  along  the  edges  of  the  map.  When  the 
mouse  cursor  is  located  on  the  map,  the  longitude  and  latitude  of  the  tip  of  the  arrow  are 
displayed  at  the  bottom  of  the  window. 

The  ground  track  map  shows  the  ground  trace  against  a  map  of  the  world.  The  outage  map 
covers  an  area  encompassing  the  ground  site  grid  entered  by  the  user  in  the  second  dialog  box 
and  plots  every  ground  site  that  was  not  covered  for  one  hundred  percent  of  the  time.  The 
symbol  used  to  plot  each  site  indicates  how  much  time  the  site  was  blacked  out.  A  different 
symbol  is  used  for  each  outage  zone,  as  indicated  on  the  map  legend.  The  spot  beam  map 
displays  spot  beams  and  possibly  loss  contours  and  the  horizon  as  described  above  on  a  map 
of  the  world.  Examples  of  all  three  types  of  maps  are  shown  in  the  appendices. 

The  user  can  expand  a  rectangular  portion  of  the  map  (“zoom  in”)  by  moving  the  mouse 
cursor  to  one  comer  of  the  desired  rectangle,  pressing  the  mouse  button,  and  dragging  a 
rectangle  to  cover  the  desired  area.  A  sample  rectangle  is  shown  in  figure  2-14.  Once  the 
zoom  operation  has  been  performed,  the  displayed  area  can  be  reenlarged  by  dragging  an 
interior  rectangle  or  shifted  slightly  by  dragging  the  rectangle  outside  of  the  map  area. 


Figure  2-14.  Map  Zooming 
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When  expanding  an  outage  zone  map,  the  user  is  given  the  opportunity  to  recalculate  the 
outage  zones  within  the  new  rectangle.  The  dialog  box  shown  in  figure  2-15  appears  on  the 
screen.  The  user  can  enter  new  step  sizes  for  the  calculation.  The  default  step  sizes  are  taken 
from  the  original  document.  Pressing  the  “Yes”  button  starts  the  recalculation.  Pressing  “No’ 
zooms  in  without  recalculating.  Pressing  “Cancel”  cancels  the  zoom  operation  altogether.  If 
the  user  presses  “Yes,”  the  recalculated  area  becomes  the  new  “zoomed  out”  state.  In  other 
words,  the  user  may  not  zoom  out  to  the  state  prior  to  the  recalculation.  (Zooming  out  is 
described  below.) 


mould  you  like  to  recalculate  with 
L3  new  step  sizes? 


No  Cancel 

>  _ 

New  latitude  step  size  (deg):  |5.0 
New  longitude  step  size  (deg):  5.0 

Note :  If  you  do  recalculate,  the  original  grid  will  no 
longer  be  available  through  the  Zoom  Out  command. 


Figure  2-15.  Outage  Zone  Recalculation  Dialog  Box 


Show  Tcxigaigmi 

This  command  is  available  if  the  active  window  is  a  map  window.  It  simply  brings  to  the 
front  of  the  screen  the  corresponding  text  window,  redisplaying  it  if  it  were  closed  or  hidden. 

Zoom  Out 

This  command  resets  the  display  area  of  the  map.  For  ground  track  and  spot  beam 
projection  documents,  this  is  the  entire  world.  For  outage  zone  documents,  this  is  the  area 
encompassing  the  ground  site  grid. 
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Show  Minute  Features 

This  command  determines  how  much  detail  will  be  displayed  on  the  map.  Selecting  this 
command  toggles  it  on  or  off;  the  on  state  is  indicated  by  the  presence  of  a  check  marie.  When 
the  feature  is  activated,  all  map  features  are  plotted,  no  matter  how  small.  When  this  feature  is 
turned  off,  the  smallest  islands,  lakes,  and  countries  are  not  displayed.  The  trade-off  is  in 
execution  time. 

2.4.6  The  Windows  Menu 

This  menu  contains  a  list  of  all  the  windows  currendy  in  memory,  including  those 
temporarily  hidden.  Selecting  a  window  from  this  menu  brings  it  to  the  front  on  the  screen, 
redisplaying  it  if  it  were  hidden. 


2.5  PRINTING  AND  OPENING  DOCUMENTS  FROM  THE  FINDER 

The  printing  of  documents  from  the  finder  is  supported  by  ESP.  To  do  this  the  user  first 
selects  the  icons  of  the  desired  files  and  then  chooses  the  Print  command  from  the  File  menu. 
The  Page  Setup  dialog  box  is  presented  first,  and  the  options  set  here  will  apply  to  all  the 
documents.  Separate  dialog  boxes  will  then  be  presented  for  each  document  for  entering  the 
final  job  settings,  e.g..  the  pages  to  be  printed  and  the  number  of  copies. 

ESP  can  also  be  launched  by  opening  any  of  its  documents  from  the  finder.  The  icons  of 
the  different  types  of  documents  are  shown  in  figure  2-1. 


2.6  ORBITAL  PARAMETERS  CONVERSION  UTILITY 

In  order  to  accommodate  different  orbital  parameter  sets,  a  conversion  utility  called  Convert 
Parameters  is  provided  on  the  diskette.  The  program's  operation  is  simple.  When  started,  the 
program  displays  the  dialog  box  shown  in  figure  2-16.  The  user  first  selects  one  of  the  five 
conversions  shown  in  the  dialog  box.  The  input  fields  in  the  dialog  box  will  change  according 
to  the  selected  conversion.  For  those  conversions  that  require  date  and  time  inputs,  the 
defaults  are  the  same  as  those  of  ESP.  After  entering  the  inputs  the  user  should  press  the 
“Convert”  button.  A  window  will  be  displayed  with  the  converted  parameter(s).  Shown  in 
figure  2-17  is  the  output  window  that  shows  the  results  corresponding  to  the  inputs  in  figure 
2-16.  The  user  can  make  as  many  conversions  as  desired  by  selecting  Convert. . .  from  the  Edit 
menu. 
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Convert... 

<8>  geosynchronous  parameters  to  orbital  elements 
O  position  and  velocity  to  orbital  elements 
O  orbital  period  to  semimajor  avis  (‘a’) 

O  mean  motion  to  ‘a’ 

O  mean  anomaly  and  ‘a’  to  time  of  perigee 
Geosynchronous  Satellite 

Subsatellite  Longitude  at  Equator:  |0.0  °  east 

Inclination:  0.0  ° 


Time  of  Equatorial  Crossing  _ 

at  Ascending  Node:  Jan  1,  1991 


00:00:00 


Figure  2-16.  Sample  Convert  Parameters  Dialog  Box 


Convert  Parameters 


Inputs 

Subsatellite  longitude  at  equator  =  0.0s  east 
Inclination  =  0.0° 

Time  of  equatorial  crossing 
at  ascending  node  =  Jan  1,  1991  00:00:00 

Orbital  Elements 

These  numbers  can  be  entered  into  ESP. 

o  -  42163.0  km 

e  -  0.0 

/  =  0.0° 

i?  =  99.41377° 

co  =0.0° 

t  =  19910101  0 

Time  of  osculating  values  =  19910101  0 


Figure  2-17.  Sample  Convert  Parameters  Output 
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SECTION  3 

ANALYTICAL  FOUNDATIONS 


3.1  OVERVIEW 

ESP  2.0  includes  perturbations  to  a  satellite’s  orbit  that  are  due  to  the  Earth’s  oblateness  (72 
effects),  atmospheric  drag,  and  solar  radiation  pressure.  Satellite  orbits  are  propagated  using  a 
classical  fourth-order  Runge-Kutta  algorithm  with  a  fixed  step  size.  This  provides  a  good 
approximation  to  a  satellite’s  position  and  velocity  vectors,  each  as  a  function  of  time.  Initial 
conditions  are  determined  from  the  values  of  the  satellite  parameters  that  are  input  by  the  user. 
Major  outputs  of  the  program  —  the  satellite  ground  track;  elevation  angle,  range  and  range 
rate  calculations;  and  spot  beam  projections  —  are  derived  from  the  satellite  position  vector  as  a 
function  of  time.  This  section  provides  the  mathematical  details  for  the  propagation  scheme  as 
well  as  all  secondary  calculations. 


3.2  FUNDAMENTALS 

To  begin,  a  review  of  some  astrodynamics  fundamentals  is  essential.  Consider  the  two- 
body  problem,  that  is,  the  simplified  scenario  in  which  the  sole  force  governing  two  masses  is 
that  of  their  mutual  gravitational  attraction.  In  this  circumstance,  the  relative  motion  is  confined 
to  a  plane  and  described  by  a  conic  section.  Furthermore,  there  are  six  integrals  or  constants  of 
motion  that  completely  define  the  behavior  of  the  particles.  These  six  constants  are  known  as 
the  orbital  elements;  they  may  assume  several  forms  and  define  the  motion  for  any  type  of 
conic.  The  classical  representation  of  these  elements  in  the  context  of  closed  orbits  (circles  and 
ellipses)  is  the  basis  for  ESP  and  subsequent  discussion. 

3.2.1  Orbital  Elements 

The  first  three  elements  are  independent  of  the  frame  of  reference  and  are  shown  in 
figure  3-1.  The  size  of  the  orbit  is  defined  by  the  semimajor  axis,  a,  which  is  one-half  the 
length  of  the  chord  passing  through  the  foci  of  the  ellipse.  The  shape  is  determined  by  the 
eccentricity,  e,  which  is  defined  as  the  ratio  of  one-half  the  distance  between  the  foci  (denoted 
by  c  in  the  figure)  over  the  semimajor  axis  (i.e.,  c/a).  For  closed  orbits,  the  eccentricity  ranges 
from  zero,  for  circles,  up  to  but  not  including  one,  for  ellipses.  The  third  element  pinpoints 
satellite  position  on  the  orbit  as  a  function  of  time;  it  is  the  time  of  closest  approach  to  the 
central  attracting  body,  known  as  the  time  of  pericenter  passage  (or  in  the  case  of  Earth- 
orbiting  satellites,  the  time  of  perigee  passage),  t. 

The  remaining  three  elements  do  in  fact  depend  on  the  reference  frame.  Typically,  when 
dealing  with  Earth -orbiting  satellites,  one  uses  the  geocentric-equatorial  coordinate  frame, 
otherwise  known  as  the  Earth- Centered  Inertial  (ECI)  coordinate  frame.  In  this  system  the 
origin  of  coordinates  is  taken  as  the  center  of  the  Earth,  with  the  x-y  plane  coincident  with  the 
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Figure  3-1.  Semimajor  Axis,  Eccentricity,  and  Time  of  Pericenter  Passage 


equatorial  plane.  The  positive  x-axis  is  in  the  direction  of  vernal  equinox1,  the  positive  z-axis 
is  along  the  axis  of  the  Earth’s  rotation  and  points  due  north,  and  they-axis  is  chosen  so  that  a 
right-handed  coordinate  system  is  formed.  It  is  important  to  note  that  this  frame  does  not  rotate 
as  the  Earth  rotates  diumally  on  its  axis.  However,  it  is  not  an  inertial  frame  in  the  true  sense 
due  to  precession2.  A  fixed  frame  of  reference  may  be  established  by  defining  the  coordinate 
system  at  some  epoch.  ESP  uses  the  mean3  equinox  and  Earth  equator  of  the  standard  epoch 
of  1950.04.  An  illustration  of  the  ECI  coordinate  frame  as  it  pertains  to  the  discussion  that 
follows  is  provided  in  figure  3-2. 


1The  vernal  equinox  is  defined  as  that  point  along  the  line  of  intersection  of  the  Earth's  equatorial  plane  and  the 
ecliptic  plane  (the  plane  of  the  Sun-Earth  motion)  at  which  the  Sun  crosses  from  south  to  north  in  its  apparent 
motion  about  the  Earth. 

2Luni-solar  precession  refers  to  the  motion  of  the  Earth's  polar  axis  about  the  surface  of  a  cone  with  period  26,000 
years.  It  is  caused  by  the  gravitational  attraction  of  the  Sun  and  Moon  on  the  Earth's  equatorial  bulge  and  results 
in  a  westward  shift  of  the  vernal  equinox  at  a  rate  of  about  50  arc-seconds  per  year.  There  is  also  an  oscillation  (or 
nodding)  of  the  Earth's  polar  axis  as  it  sweeps  out  this  cone.  This  effect  is  imposed  by  the  gravitational  attraction 
of  the  Moon,  has  period  18.6  years,  and  is  referred  to  as  nutation.  In  addition,  planetary  precession,  caused  by  the 
gravitational  attraction  of  the  planets,  results  in  a  slow  shift  in  the  orientation  of  the  ecliptic  plane. 

'Mean  equinox  and  equator  neglects  the  relatively  short-period  variations  caused  by  nutation. 

4The  standard  epoch  of  1950.0  is  defined  as  the  beginning  of  the  Besselian  year,  or  the  Julian  ephemeris  date 
2433282.423357.  The  reference  planes  are  fixed  at  this  instant.  The  time  reference  epoch  that  appears  in  the 
calculations  of  later  sections  (sections  3.3.3  and  3.6)  is  the  nearest  even  date,  which  is  zero  hours  on  1  January 
1950,  or  the  Julian  ephemeris  date  2433282.5. 
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Figure  3-2.  Euler  Angles 


Now  that  the  frame  of  reference  has  been  established,  the  remaining  three  elements  may  be 
defined.  First  note  that  the  equatorial  plane  and  the  orbital  plane  (the  plane  of  the  orbiting 
body's  motion)  intersect  in  a  line  (when  the  orbital  plane  is  different  from  the  equatorial  plane), 
which  is  referred  to  as  the  line  of  nodes.  The  ascending  (descending)  node  is  that  point  along 
tire  line  of  nodes  where  the  satellite  crosses  from  south  (north)  to  north  (south)  in  its  orbit 
The  angle  in  the  equatorial  plane  between  the  positive  x-axis  and  the  direction  of  the  ascending 
node  is  the  longitude  of  the  ascending  node,  Ll.  (For  the  case  in  which  the  orbital  plane  and 
tiie  equatorial  plane  are  the  same,  Q  is  arbitrary.)  The  angle  in  the  orbital  plane  measured  from 
the  direction  of  the  ascending  node  to  the  direction  of  perigee  is  the  argument  of  perigee,  (O. 
Finally,  the  angle  between  the  equatorial  plane  and  the  orbital  plane  is  the  inclinations,  i.  These 
three  angles  collectively  are  known  as  the  Euler  angles. 

3.2.2  Osculating  Orbit 

If  the  basic  premise  of  two-body  motion  held  in  real  life,  then  it  would  be  a  simple  task  to 
determine  the  satellite  position  and  velocity  vectors  as  functions  of  time.  (This  would  be 
accomplished  exactly  according  to  the  method  provided  in  section  3.5.)  Unfortunately,  there  is  a 
complex  interplay  of  forces  that  affects  the  motion  of  any  body  in  space.  In  actuality,  the  six 
elements  described  above  are  not  constants  but  functions  of  time,  and  the  true  motion  of  a  body 
orbiting  the  Earth  is  not  a  perfect  ellipse  but  some  irregularly  shaped  path.  At  any  instant  the  true 
orbit  is  tangent  to  a  perfectly  elliptical  orbit,  known  as  the  osculating  orbit,  with  corresponding 
orbital  elements  identical  at  that  instant  The  instant  is  thus  referred  to  as  the  time  of  the 


5The  inclination  is  often  equivalently  defined  as  the  angle  between  the  positive  2-axis  and  the  orbit  normal  vector 
(the  vector  n  in  figure  3-2). 
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osculating  values  (or  frequently,  the  epoch  time).  The  tangential  orbit  may  be  regarded  as  that 
orbit  which  would  result  if  all  perturbations  were  turned  off  at  that  instant.  Figure  3-3  depicts  the 
relationship  between  the  actual  and  osculating  orbits. 


Given  the  six  orbital  elements  of  the  osculating  orbit  at  some  particular  instant,  the  position 
and  velocity  are  relatively  easily  determined  at  that  instant.  Subsequent  positions  and  velocities 
can  then  be  estimated  using  these  initial  conditions  and  the  equations  of  motion  in  conjunction 
with  an  integration  scheme. 

3.2.3  Additional  Parameters 

In  addition  to  the  orbital  elements  discussed  above,  there  are  six  other  parameters  of  direct 
interest.  The  first  three  are  illustrated  in  figure  3-4.  The  true  anomaly,/,  is  the  angle  in  the 
orbital  plane  measured  from  the  direction  of  perigee  (often  referred  to  as  the  line  of  apsides  or 
apsidal  line)  to  the  satellite  position  vector,  r.  The  eccentric  anomaly,  E,  may  be  obtained  from 
the  true  anomaly  geometrically  as  follows.  First,  an  auxiliary  circle  with  radius  equal  to  the 
semimajor  axis  a  is  circumscribed  about  the  ellipse.  Next,  a  perpendicular  is  dropped  from  the 
satellite  position  vector  to  the  major  axis  and  extended  to  the  point  P  on  the  circle.  The 
eccentric  anomaly  is  the  resultant  angle  formed  between  the  line  of  apsides  and  the  line  segment 
OP.  There  are  also  algebraic  relationships  between  /  and  E,  however  these  are  not  of  concern 
in  the  present  analysis.  The  third  parameter  shown  in  figure  3-4  is  the  semiminor  axis,  b;  it  is 
one-half  the  length  of  the  chord  passing  through  the  center  of  the  ellipse  and  perpendicular  to 
the  major  axis.  The  semiminor  axis  is  related  algebraically  to  the  semimajor  axis  and 
eccentricity  as  follows: 


(3-1) 
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Figure  3-4.  Additional  Parameters 


The  remaining  three  parameters  needed  are  defined  algebraically.  The  time  required  for  a 
satellite  to  make  one  complete  pass  through  its  orbit,  or  the  orbit  period,  P,  is  a  function  of  the 
semimajor  axis  a  and  the  gravitational  parameter  p: 


For  an  Earth-orbiting  satellite. 


(3-2) 


p  -  G{me  +  ms) 

where  G  is  the  universal  gravitational  constant,  and  me  and  ms  are  the  masses  of  the  Earth  and 
satellite,  respectively.  Since  the  mass  of  the  satellite  is  negligible  in  comparison  with  the  mass 
of  the  Earth, 

p  ~  Gme  =  398600.45  km3  /  sec2 

Note  that  equation  3-2  can  be  rearranged  to  determine  the  semimajor  axis  associated  with  a 
given  orbit  period 


l 


r 


The  average  angular  motion  (radians  per  unit  time)  or  mean  motion ,  n,  is  defined  in  terms 
of  the  orbit  period  by 


The  final  parameter  is  the  mean  anomaly,  M.  Regardless  of  the  type  of  motion,  M  is 
always  a  function  of  time.  It  is  defined  in  terms  of  the  mean  motion  n  and  the  time  of 
pericenter  passage  x.  Explicitly, 


M(t)  =  n(t-  t) 

M(t)  may  be  interpreted  geometrically  as  the  angle  measured  from  the  line  of  apsides  to  the 
position  vector  obtained  at  time  t  if  the  satellite  were  moving  with  constant  angular  velocity  n. 
Often  the  time  of  the  osculating  values,  tosc  and  the  time  of  pericenter  passage,  x,  are  provided 
implicitly  by  the  single  parameter  M(tosc). 


3.3  EQUATIONS  OF  MOTION 


In  the  case  of  two-body  motion,  the  sole  force  exerted  on  the  two  mass  particles,  mi  and 
m2,  is  that  of  their  mutual  gravitational  attraction.  In  accordance  with  Newton’s  Law  of 
Universal  Gravitation,  the  magnitude  of  this  force  is  directly  proportional  to  the  product  of  the 
masses  and  inversely  proportional  to  the  square  of  the  distance  between  them.  The  direction  of 
the  force  exerted  by  m2  on  mi  is  along  the  line  of  sight  from  mi  to  m2  (and  oppositely  oriented 
is  the  force  exerted  by  mi  on  m2).  Based  on  Newton’s  Second  Law6,  the  motion  of  each  of 
the  masses  with  respect  to  a  common  origin  O  (refer  to  figure  3-5)  is  expressed  algebraically  as 
follows: 


..  _Gmxm2  ni+ir 
iTi - -5—  (“I)  - 


1  =  1,2 


where  r,-  represents  the  position  of  mass  mi  relative  to  origin  O,  r  is  the  position  of  m2  relative 
to  mi  (i.e.,  r  =  r2  -  ri),  and  r  is  the  magnitude  of  the  vector  r. 

Subtraction  of  the  case  i  =  2  from  the  case  i  =  1  yields  the  following  equation  that  completely 
describes  the  relative  motion  of  the  particles: 


(3-4) 


6 Newton's  Second  Law  states  that  the  rate  of  change  of  momentum  is  proportional  to  the  force  impressed  and  is  in 
the  same  direction  as  that  force  (lj.  For  a  body  with  constant  mass  m,  it  is  expressed  mathematically  as 
F  =  na,  where  F  is  the  net  force  exerted  on  the  body  and  a  is  the  acceleration  of  the  body. 
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Figure  3-5.  Two-body  Geometry 


This  second-order  equation  is  often  expressed  equivalently  as  a  system  of  two  first-order 
equations: 


dx__  dv  ^  _ 

dt  dt  ~r$ 


where  v  is  the  velocity  of  m2  relative  to  m\.  In  the  presence  of  additional  forces,  equation  3-4 
becomes: 


r  = 


(3-5) 


where  F  represents  the  net  disturbing  force  acting  on  the  system.  F  is  simply  the  sum  of  all 
the  individual  disturbing  forces,  i.e.. 


F  =  fj  -t-f2  +  f3+... 


where  each  fj  is  an  individual  disturbing  force.  Since  these  individual  disturbing  forces  cause 
the  motion  to  deviate  from  a  perfect  conic,  they  are  commonly  referred  to  as  perturbations. 

ESP  2.0  takes  into  account  the  effects  of  three  perturbations,  namely,  the  oblateness  of  the 
Earth  (J2  effects),  atmospheric  drag,  and  solar  radiation  pressure.  The  mathematical  model  for 
each  is  presented  in  detail  in  the  subsections  that  follow.  As  a  notational  convenience,  they 
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will  be  denoted  by  fi,  f2,  and  f3,  respectively.  Also,  the  satellite  position  and  velocity  vectors 
with  respect  to  the  ECI  coordinate  frame  will  be  denoted  by 


respectively. 

3.3.1  Perturbations  Due  to  J2  Effects 


The  major  force  influencing  the  motion  of  an  Earth-orbiting  satellite  is  the  Earth’s 
gravitational  field  If  the  Earth  were  a  perfectly  smooth  sphere  with  a  uniform  mass 
distribution,  then  equation  3-4  would  completely  describe  the  motion  resulting  from  this  force. 
In  actuality,  the  Earth  is  pear-shaped,  meaning  it  is  more  massive  below  the  equatorial  belt  than 
above.  Under  the  assumption  that  the  Earth  is  axially  symmetric,  the  additional  force  exerted 
on  a  satellite  due  to  the  Earth’s  uneven  mass  distribution  is  determined  from  the  Earth’s 
gravitational  potential  function 


where  J/c  =  empirically  determined  constants 

Req  -  equatorial  radius  of  the  Earth  =  6378.14  km 
Pk  =  Legendre  polynomials,  k-  0,1, 2, ... 

In  the  above  expression,  the  first  two  terms  are  dominant,  a  result  of  the  fact  that 

(a)  the  coefficients  are  all  very  close  to  zero,  and  become  significantly  more  so 
as  k  increases  (as  an  indication,  Ji  is  smaller  than  by  three  orders  of 
magnitude) 

(w 


and 


(c) 


<1 


Neglecting  all  subsequent  terms  results  in 


V(r,i)= 
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where 


72  =  zonal  harmonic  of  degree  two  =  0.00108263  [3] 


and  /i(u)  =  i(3u2-l) 

The  force  per  unit  mass  acting  upon  the  satellite  due  to  72  is  found  by  taking  the  gradient  of 

V*  =  V  -  — 
r 

with  respect  to  the  position  vector  r,  i.e., 


'dV *' 

T 

' dV * 
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dx 
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dz 
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Note  that  the  contribution  to  the  net  force  that  is  a  result  of  the  first  term  of  V,  /Mr,  is  included 
as  the  first  term  in  equation  3-5. 

3.3.2  Perturbations  Due  to  Atmospheric  Drag 

Atmospheric  drag  is  an  important  force  influencing  the  motion  of  satellites  in  low-Earth 
orbit.  The  general  effect  of  this  force  is  that  it  causes  energy  to  dissipate.  Drag  is  directly 
proportional  to  atmospheric  density  which  in  turn  decreases  with  altitude  above  the  Earth’s 
surface.  It  is  also  a  function  of  a  satellite’s  cross-sectional  area-to-mass  ratio.  For  any  given 
system,  the  lower  the  orbit  the  greater  the  impact  of  drag.  In  general,  orbits  between  100  km 
and  200  km  above  the  Earth  will  decay  in  a  matter  of  a  few  revolutions,  whereas  orbits  higher 
than  500  km  will  have  lifetimes  measured  in  years.  The  impact  of  drag  is  essentially  negligible 
for  orbits  that  are  higher  than  1000  km  above  the  Earth. 
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Mathematically,  the  force  exerted  on  a  satellite  due  to  drag  is  expressed  as  follows: 


where  Co 

Ad 

M 

P 

Vrel 

vrel 


.  1  CnA[) 

coefficient  of  drag,  dependent  on  vehicle  shape  and  attitude;  several 
sources  recommend  a  value  of  2.0 

effective  satellite  area  for  drag  (cross-sectional  area  perpendicular  to 
the  direction  of  motion) 
satellite  mass 

atmospheric  density  corresponding  to  satellite  altitude 
velocity  of  the  satellite  relative  to  the  rotating  atmosphere 
magnitude  of  \re[ 


ESP  uses  a  static  atmospheric  density  model  based  on  the  1976  U.  S.  Standard 
Atmosphere  and  adapted  from  the  JPL  tool  ASAP7.  This  tabular  model  provides  atmospheric 
density  as  a  function  of  altitude,  for  altitudes  ranging  from  86  km  to  1000  km.  Altitudes  lower 
than  86  km,  though  highly  unlikely,  are  assumed  to  have  density  equivalent  to  that  at  86  km. 
Altitudes  higher  than  1000  km  are  assumed  to  have  zero  density.  The  densities  of  altitudes 
between  values  in  the  table  are  determined  by  linear  interpolation. 

The  altitude,  K  of  the  satellite  at  any  instant  is  given  by 


h~r-Re 


where  r  is  the  magnitude  of  the  position  vector,  and  Re  is  the  radius  of  the  Earth.  ESP  models 
the  Earth  as  an  oblate  spheroid  in  the  determination  of  Re.  In  general,  the  radius,  £,  of  an 
oblate  spheroid  is  given  by 


where  ae  and  et  are  the  semimajor  axis  and  eccentricity,  respectively,  of  an  elliptical  cross- 
section  through  the  poles,  and  o  is  the  geocentric  latitude.  An  expression  for  the  eccentricity  of 


7The  Artificial  Satellite  Analysis  Program  (ASAP)  [IS]  is  a  sophisticated  trajectory  propagator  useful  for 
analyzing  orbit  behavior  over  a  period  of  several  months.  It  includes  perturbations  on  the  spacecraft  orbit  due  to 
non-sphericity  (up  to  a  40  x  40  field)  of  the  central  body,  luni-solar  effects,  drag  and  solar  radiation  pressure,  and 
uses  an  eighth-order  Runge-Kutta  integration  scheme  with  variable  step  size  control.  The  software  runs  on  a  PC 
and  takes  as  input  a  lengthy  column-formatted  ftle.  Major  outputs  are  time  histories  of  the  orbital 
elements/satellite  state  in  a  variety  of  representations. 
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this  cross-section  in  terms  of  the  semimajor  and  semiminor  axes  is  found  by  rearranging 
equation  3-1: 


For  the  Earth,  ae  and  be  are  the  equatorial  and  polar  radii,  respectively,  as  illustrated  in  figure 
3-6.  Their  values  are  given  by 

ae  =  6378  km 
be  *  6357  km 


Figure  3-6.  Earth  Cross-section 


which  leads  to  the  following  widely  referenced  ([2],  [8],  [9],  [10],  [1 1  j,  [15])  value  for  ee: 

ee  =  0.08182 

(Note  that  this  result  is  obtained  using  more  precise  values  for  ae  and  be  than  the 
approximations  indicated  above.)  Also,  at  any  particular  time  the  latitude  5  corresponding  to 
the  satellite  position  vector  r  is  given  by 

S  =  arcsin 

as  shown  in  section  3.6. 
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Finally,  the  velocity  of  the  satellite  relative  to  the  rotating  atmosphere  is  needed.  This  is 
obtained  by  direct  application  of  the  Coriolis  Theorem.  This  theorem  provides  a  relationship 
between  a  vector  in  a  fixed  coordinate  system  F  and  the  same  vector  in  another  coordinate 
system  R,  rotating  with  respect  to  F  with  angular  velocity  vector  co.  For  any  vector  a  in  F,  this 
relationship  is  expressed  as  follows: 


where  in  general, 


da 

+(oxa 

dtR 


denotes  the  time  derivative  of  a  with  respect  to  the  coordinate  frame  A. 


Now  the  Earth  is  rotating  in  the  ECI  frame  with  mean  angular  velocity 


“e  = 


0 

0 

L®*. 


where 


Q)e  =15°.041067178  per  hour 

Application  of  the  above  theorem  to  the  position  vector  r  in  the  fixed  or  ECI  coordinate  frame 
yields  the  velocity  of  the  satellite  relative  to  the  rotating  atmosphere,  as  follows: 


-G)exr 


=  v  -  (oe  x  r 


Vjc  +yo>e 
Vy-XO)e 
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3.3.3  Perturbations  Due  to  Solar  Radiation  Pressure 


Solar  radiation  pressure  is  caused  by  the  light  energy,  or  photons,  of  the  Sun.  Depending 
(hi  where  the  satellite  is  in  its  orbit,  it  may  add/subtract  energy  to/from  the  orbit  As  one  might 
expect,  the  force  exerted  on  a  satellite  due  to  solar  radiation  pressure  is  a  function  of  the 
position  of  the  Sun  relative  to  the  Earth.  Like  drag,  it  is  also  a  function  of  a  satellite’s  cross- 
sectional  area-to-mass  ratio  (although  the  effective  satellite  area  for  solar  radiation  may  very 
well  differ  from  that  for  atmospheric  drag).  In  extreme  cases  where  the  area-to-mass  ratio  is 
very  large  (e.g.,  balloon  satellites)  the  effect  of  solar  radiation  is  to  “blow”  the  satellite  off 
course,  a  phenomenon  referred  to  as  “solar  sailing”  and  one  that  may  be  desirable  for 
spacecraft  navigation  or  attitude  control  This  force  is  modelled  mathematically  as  follows: 


*3  =  - 


YPsAs_  y 
M  s 


=  reflectivity/absorption  constant,  0  £  y<,  2 
<  1  for  translucent  material 

=  1  for  perfectly  absorbent  material 

=  2  for  perfectly  reflective  material 

=  radiation  pressure  (force  per  unit  area)  on  a  perfectly  absorbing 
surface  at  1  astronomical  unit8 
=  4.4  x  lO-3  kg/km-sec2 

=  effective  satellite  area  for  solar  radiation  pressure  (cross-sectional 
area  exposed  to  the  Sun) 

=  satellite  mass 

=  unit  vector  to  the  Sun  in  the  ECI  frame,  i.e.,  direction  of  the  Sun 
with  respect  to  the  Earth 

The  direction  of  solar  radiation  pressure  is  actually  assumed  to  be  along  the  line-of-sight 
from  the  Sun  to  the  satellite,  or  r  -  rs.  Since  the  Sun  is  located  at  a  large  distance  from  both  the 
Earth  and  a  satellite  orbiting  the  Earth,  r  -  rs  *»  -  rs.  The  position  vector  to  the  Sun  in  the  ECI 
frame,  rs,  is  obtained  by  considering  the  Sun  as  a  body  orbiting  the  Earth.  Given  the  orbital 
elements  of  the  Sun  at  time  t,  rs (f)  is  determined  exactly  according  to  the  method  outlined  in 
section  3.5.  Reference  21  provides  the  following  relationships  for  the  mean  orbital  elements  of 
the  Sun  with  respect  to  the  mean  equinox  and  Earth  equator  of  1950.0  (distance  in  kilometers, 
and  angles  in  degrees): 

as  =  1.49597927  x10s 

es  =  1.67301085  xl0~2  -4. 1926  xl0~5r- 1.26  xlO'7T2 


where  y 

P s 

As 

M 

A 

*5 


8 An  astronomical  unit  (A.  U.)  is  defined  as  the  mean  distance  from  the  Earth  to  the  Sun,  or  approximately  ISO 
million  kilometers. 


37 


is  =  23.4457888616  - 1.30141669  x  lO~2T-9.44SxlO~1  T2  +  5.000  x  10~7r3 
Q)s  =  282.08053  +  3. 2328  x  10"1  T + 1. 5  x  10^  T2 

Ms  =  358.000682  +  9.856002623  x  10-1d  - 1.550000 x  10-4 T2  - 3.3333 x  lO^T3 

where  d  =  number  of  ephemeris  days  from  reference  epoch  1  January  1950  to  time  t 

and  T  =  number  of  Julian  centuries  from  reference  epoch  1  January  1950  to  time  t 
=  <036525 

Note  that  the  longitude  of  the  ascending  node,  Qs,  is  zero  since  the  positive  Jt-axis  is  defined  as 
the  direction  of  the  vernal  equinox.  Also  note  that  the  first  four  elements  listed  above  are 
essentially  fixed  for  the  duration  of  the  run.  Thus,  they  are  calculated  once  at  starting  time.  The 
mean  anomaly,  however,  determines  the  location  of  the  Sun  in  its  orbit  relative  to  the  Earth  as  a 
function  of  time.  Therefore,  it  must  be  recalculated  with  each  evaluation  of  the  function  f3. 


In  the  case  in  which  the  Earth  blocks  or  “shadows”  the  Sun  from  the  satellite,  there  is  no 
solar  radiation  pressure  effect.  Shadowing  may  occur  if  the  angle  6  subtending  the  vectors  r 
and  rs  is  greater  than  90°  (refer  to  figure  3-7).  0is  calculated  from  the  dot  product  as  follows: 


0  =  arccos 


If  0  does  in  fact  exceed  90°,  the  Earth  will  block  the  line-of-sight  vector  from  the  Sun 
to  the  satellite  when 

rsin0<  Rfotn 

where  Ratm  denotes  the  radius  of  die  Earth  plus  the  obscuring  atmosphere  and  is  approximated 
by  Re  +  90  km  [15]. 


Figure  3-7.  Test  for  Shadowing 
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3.4  INTEGRATION  SCHEME 

The  equations  of  motion  in  section  3.3  have  no  closed-form  solution.  However, 
approximations  to  a  satellite’s  position  and  velocity  vectors  as  functions  of  time  can  be 
obtained  using  numerical  methods.  The  scheme  used  by  ESP  is  a  classical  fourth-order 
Runge-Kutta  algorithm. 

To  develop  a  general  understanding  of  the  theory  behind  this  algorithm,  consider  the 
equations  for  two-body  motion 


dr 

dt 


v 


(3-7) 


with  initial  conditions  at  time  t0,  r(t0)  and  v(r0).  The  position  and  velocity  at  some  later  time 
tQ  +  At  can  be  found  by  expanding  r (t)  and  v(r)  each  as  a  Taylor  series  about  the  value  t0.  As 
an  example,  first-order  approximations  are  given  by 


r(r0  +  At)  =  r(t0  )  +  ^|  At  +  O(At) 
dt 


v(r0+AO  =  v(r0)  + 


At+0(At)2 


where  terms  of  the  order  (At)2  and  higher  are  neglected.  (Note  that  this  example  serves  for 
illustration  purposes  only.  As  section  3.2  indicates,  the  six  orbital  elements  explicitly  define 
two-body  motion;  as  such,  an  integration  scheme  is  neither  needed  nor  desirable.)  In  general, 
an  n*  older  algorithm  neglects  terms  of  order  (At/1  * 1  and  higher.  Accuracy  is  improved  by 
using  a  higher  order  algorithm  or  a  smaller  time  increment,  or  both.  However,  when  moving 
to  higher  orders,  the  derivatives  on  the  right-hand  side  of  the  above  equations  can  become  quite 
complex.  Even  with  a  second-order  algorithm,  the  third  term  in  the  expansion  for  the  velocity 
vector  includes 


d2\  _  dg  dr 

dt 2  ,  dr  dt . 

t0 

which  involves  calculating  a  3  x  3  matrix.  Moving  to  a  third-order  algorithm  would  involve  a 
3-dimensional  array,  and  the  computations  become  increasingly  more  complicated  with  orders 
beyond.  These  cumbersome  calculations  can  be  avoided  by  considering  Taylor  series 
expansions  of  appropriate  vector  sums  and  making  substitutions.  The  coefficients  in  these 
sums  can  be  chosen  in  such  a  way  that  the  higher  order  algorithms  all  follow  the  same  general 
format  An  underlying  goal  in  making  these  substitutions  is  to  minimize  the  number  of 
function  evaluations  (i.e.,  evaluations  of  the  vector  g),  otherwise  known  as  the  number  of 
stages.  Inspection  of  the  equations  of  motion  with  just  three  perturbations  included  illustrates 
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the  desirability  of  attaining  this  goal.  The  details  for  these  conversions  are  themselves  long  and 
tedious  and  are  deferred  to  reference  3. 

The  general  form  of  equations  3-7  does  not  hold  when  including  atmospheric  drag  and 
solar  radiation  pressure  effects.  In  particular,  the  vector  g  is  no  longer  strictly  a  function  of  r 
but  now  also  a  function  of  v  and  t.  This  problem  is  handled  by  defining  the  following  6x1 
state  vector 


so  that 


ds  v 
-  =  fM 

Given  the  initial  value  s0  defined  at  time  t0,  the  classical  four-stage,  fourth-order  algorithm  for 
approximating  the  state  at  time  t0  +  At  is  as  follows: 

s  =  sc  + — (k0  +  2kj  +  2k2  +  k3)  +  O(At)5 
6 


where 


k0  =  f(t0,s0) 


*1 


k2 


^3-t(t0  +  At,s0  +  AA2) 

The  new  value  for  s  corresponding  to  time  t0  +  At  then  serves  as  the  initial  condition  for 
obtaining  the  state  at  time  t0  +  2 At,  and  the  process  repeats. 

ESP  2.0  uses  a  fixed  step  size  (At)  of  approximately  one  minute.  The  error  introduced  due 
to  a  fixed  step  size  and  truncation  of  the  series  expansions  after  five  terms  was  found  to  be 
minimal,  as  addressed  in  section  4.4. 
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3.5  INITIAL  CONDITIONS 


This  section  provides  the  details  for  determining  satellite  position  and  velocity  vectors  at 
some  instant  given  the  six  orbital  elements  at  that  instant  ( a ,  e,  i,  A  co,  t,  and  tosc,  or 
equivalently,  a,  e,  i,  A  to,  and  M(tosc) ).  The  following  method  is  used  to  obtain  initial  values 
for  the  integration  scheme  as  well  as  for  updating  the  position  vector  to  the  Sun  as  a  function  of 
time. 

First  consider  the  (x',  y\  z')  coordinate  system  obtained  by  rotating  the  ECI  coordinate 
frame  through  the  three  Euler  angles,  A  o>,  and  i  (see  figure  3-8).  In  the  rotated  system,  the 
x’-y’  plane  corresponds  to  the  orbital  plane  with  the  positive  x'  direction  toward  perigee.  The 
position  vector  in  the  rotated  frame  may  be  expressed  in  terms  of  the  eccentric  anomaly, 
semimajor,  and  semiminor  axes,  as  follows: 


x 


Figure  3-8.  Rotated  Coordinate  System 


The  velocity  vector  is  calculated  by  differentiating  r  with  respect  to  time,  which  results  in 


V  =  ' 


na 


-a  sin  E1 
bcosE 
0 


since 


dE  _  na 
~dt~  r 

The  z'  component  of  these  vectors  is  naturally  zero  since  both  lie  in  the  plane  of  motion. 

The  eccentric  anomaly  is  related  to  the  mean  anomaly  and  eccentricity  by  Kepler’s 
equation 


E  =  M  +  esin  E 

which  is  transcendental  in  E.  It  is  easily  shown  that  this  equation  has  a  unique  solution,  and  a 
variety  of  numerical  solution  techniques  exist.  Reference  3  provides  a  simple  proof  as  well  as 
several  methods  that  may  be  used  to  solve  for  E.  ESP  uses  Newton’s  method,  which  is  an 
iterative  technique.  An  initial  guess  is  taken  as  M  +  e  sinAf;  subsequent  values  are  determined 
from  the  relation 


Ek+l  ~ 


Mk) 

f'(Ek) 


(3-8) 


where 


f(E)  =  E-es in  E-M 


until  successive  values  agree  to  some  desired  level  of  accuracy.  (Accuracy  is  further  addressed 
in  section  4.)  This  technique  actually  works  quite  well,  often  converging  within  four 
iterations.  An  exception  occurs  when  e  is  very  close  to  one  and  E  is  close  to  zero  or  lit.  In 
this  case,  the  denominator  in  equation  3-8  approaches  zero  and  convergence  difficulties  may  be 
encountered.  However,  with  an  initial  guess  taken  as  n,  this  problem  should  be  avoided  [2]. 

The  position  and  velocity  in  the  ECI  frame  are  then  obtained  by  left  multiplying  the  position 
and  velocity  in  the  rotated  frame  by  the  rotation  matrix 


R  = 


cos  £2  cos  a)  -  sin  £2  sin  co  cos  i 
sin  jRcos  co  +  cos  £2  sin  co  cos/ 
sin  (0  sin  i 


-  cos  £2  sin  co  -  sin  £2  cos  co  cos  i 

-  sin  Q  sin  co  +  cos  £2  cos  co  cos  / 

cos  co  sin/ 


sin  £2  sin  i 
-cos  £2  sin  i 
co  si 
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whose  entries  are  the  direction  cosines  of  the  nine  angles  formed  between  the  axes  of  the  two 
coordinate  systems. 


3.6  DETERMINATION  OF  SATELLITE  GROUND  TRACK 

At  any  particular  instant  the  satellite  position  vector  intersects  the  surface  of  the  Earth  at  a 
point  referred  to  as  the  subsatellite  point.  The  satellite  ground  track  is  generated  by  connecting 
a  discrete  set  of  subsatellite  points  that  result  as  the  satellite  moves  in  its  orbit  over  time.  The 
location  of  a  particular  subsatellite  point  in  the  ECI  frame  is  obtained  by  projecting  the  satellite 
position  onto  the  Earth  and  can  be  described  by  the  two  angles,  asat  and  shown  in  figure 
3-9.  These  angles  are  commonly  referred  to  as  the  right  ascension  and  declination ,  respectively 
They  are  related  to  the  satellite’s  Cartesian  coordinates  as  follows: 


Figure  3-9.  Subsatellite  Right  Ascension  and  Declination 


Note  that  the  geocentric  latitude  is  synonymous  with  the  declination  8sat.  However,  the 
longitude  is  measured  from  the  prime  meridian  (the  0*  meridian  passing  through  Greenwich, 


43 


England).  The  longitude,  Ysat .  of  the  subsatellite  point  can  be  obtained  from  the  subsatellite 
and  prime  meridian  right  ascensions,  (Xsat  and  PM,  respectively,  by  the  following  relation 

Ysat  ~  asat  ~  PM 


Reference  9  provides  the  following  approximation  for  the  ephemeris  position  of  PM  at  time  t  in 
the  standard  coordinate  system  of  1950.0: 


PM  =  99 
“.87  +  24 (Ok)(d) 


where  as  before, 

(oe  =  mean  rotation  rate  of  the  Earth  =  15°.041067178  per  hour 
and  d  =  number  of  ephemeris  days  from  reference  epoch  1  January  1950  to  time  t 

Also  note  that  the  typical  problems  that  occur  with  standard  arctangent  functions  (e.g.  when  the 
argument  is  equal  to  90“  or  270”  or  falls  in  quadrants  II  or  HI)  are  avoided  by  using  a  “smart” 
arctangent  function,  as  described  in  section  4.6. 


3.7  RANGE  AND  RANGE  RATE  CALCULATIONS 

Range  is  defined  as  the  distance  from  a  given  ground  site  location  to  the  satellite.  At  any 
particular  instant,  this  distance  is  determined  as  a  function  of  the  position  vectors  in  the  ECI  frame 
to  the  site  and  the  satellite.  The  position  vector  to  the  site  is  defined  in  terms  of  its  latitude  and 
longitude  coordinates.  As  described  above,  the  declination,  5site}  is  equivalent  to  the  latitude,  and 
the  right  ascension,  aJlW,  is  a  function  of  the  longitude  and  the  prime  meridian  location: 

S Site  =  site  latitude 
asite  =  PM  +  site  longitude 

Also,  the  magnitude  of  this  vector,  rsue,  is  identical  to  the  radius  of  the  Earth  at  the  given 
location;  it  is  determined  as  a  function  of  the  site  latitude  using  equation  3-6.  The  vector  to  the 
site  in  the  ECI  frame  is  then  obtained  using  spherical  coordinates. 


xsite 

'cos(6site)cos(asite)' 

Tsite  ~ 

Ysite 

=  rsite 

cos(Ssite)sin(asiJe) 

Mite . 

- 1 

/-—N. 

£ 

£ 

V) 

_ 1 

As  shown  in  figure  3-10,  the  line-of-sight  vector  from  the  ground  site  to  the  satellite  is 
calculated  as  the  vector  subtraction, 

d  =  r  -  rsite 
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Figure  3-10.  Range  Geometry 
The  range  is  then  simply  the  magnitude  of  the  vector  d,  i.e., 

1 

M  =  [(x-*s tie)2  +C y~  Vsite )2  +(*  -  Zn7<)2]2 

Range  rate  is  calculated  as  the  time  derivative  of  the  range,  or 

d[|dQ  dd 

dt  ~  |d| 

with  the  time  derivative  of  the  vector  d  given  by 


x  xsite 

vx  ■*"  rsite  co^Ssite)s\n(tXsile)(Oe 

d  = 

y-ysite 

II 

vy  ~  r^te  cos( <5^ )  cos(  )ae 

j  ~  *site  . 

i — 

< 

i 

o 

1 _ 

since 

Aasite  =  (OeAt  =>  =  6)e 

AS  site  ~  0  =>  =  0 
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3.8  VISIBILITY  TEST 


In  general,  visibility  between  a  particular  ground  site  location  and  a  single  satellite  is 
possible  when  the  line-of-sight  connecting  them  is  above  the  horizon.  In  other  words,  the 
elevation  angle  measured  from  the  horizon  to  the  vector  d  must  be  greater  than  or  equal  to  zero. 
Figure  3-1 1  illustrates  the  geometry.  The  vector  a  is  normal  to  the  horizon  and  intersects  the 
x-y  plane  in  the  geodetic  latitude,  L.  L  is  related  to  the  site  declination  (otherwise  known  as  the 
geocentric  latitude),  by  the  following  expression: 


L 


=  arctan 


Figure  3-11.  Elevation  Angle  Geometry 

where  ee  is  the  eccentricity  of  an  elliptical  cross-section  through  the  Earth  (equal  to  0.08182)  as 
defined  in  section  3.3.2.  (For  a  derivation,  see  reference  10.)  Thus,  the  planar  triangle  OPQ 
has  known  angles  |5jj{<|,  \L\  -  |5si-(e|,  and  180"-|L|.  The  vector  b  is  therefore  determined  as 
follows: 

’cos(aJjfe)’ 
b  =  b  sin (ctsug) 

0 
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where  o&te  is  the  site  right  ascension  and  the  magnitude,  b,  is  calculated  from  the  planar 
triangle  law  of  sines  as 


b  _  tag  sinfllj-fotJ) 

sin(180*-jL|) 


(3-9) 


Note  that  when  |L|  =  0  then  b  =  0;  thus,  the  singularity  in  equation  3-9  is  easily  avoided.  The 
vector  a  is  then  given  by  the  vector  subtraction 

a  =  b  “  *site 

Finally,  the  angle  yf  between  the  two  vectors  a  and  d  is  calculated  from  the  dot  product  as 

yf  =  arccos 


The  elevation  angle  is  simply  yf- 90°.  A  negative  value  indicates  that  the  satellite  is  below  the 
horizon. 

In  the  outage  portion  of  the  program,  the  user  inputs  a  minimum  elevation  angle  below 
which  satellites  and  ground  sites  cannot  see  each  other.  Coverage  is  obtained  at  a  given  ground 
site  location  only  when  the  elevation  angle  measured  from  this  site  is  above  this  minimum 
elevation  angle  for  the  user-specified  minimum  number  of  satellites  in  the  constellation. 


3.9  SPOT  BEAM  PROJECTIONS 

The  term  spot  beam  refers  to  a  satellite  beam  that  is  narrow  in  comparison  to  an  Earth 
coverage  beam.  A  spot  beam  projection  is  defined  as  the  contour  of  intersection  of  such  a  beam 
and  the  surface  of  the  Earth.  ESP  assumes  that  the  beam  is  a  circular  cone.  The  user  specifies  a 
pointing  location  for  the  beam  centroid,  i.e.,  longitude  and  latitude  coordinates  on  the  surface  of 
the  Earth,  and  a  beam  width,  defined  as  the  angle  that  subtends  the  full  cone.  The  technique 
used  for  generating  these  projections  is  based  on  the  one  described  in  reference  1  with 
modifications  which  are  described  below.  The  basic  geometry  is  illustrated  in  figure 
3-12. 


The  idea  behind  the  technique  is  to  determine  a  sufficient  number  of  discrete  points  along 
the  intersection  of  the  cone  and  the  Earth’s  surface  such  that,  when  these  points  are  connected, 
the  resultant  curve  will  appear  smooth.  In  the  terminology  of  reference  1,  these  intersection 
points  are  referred  to  as  generator  points.  The  coordinates  of  an  individual  generator  point  are 
derived  from  the  angle  de  on  the  surface  of  the  Earth  (which  may  alternatively  be  thought  of  as 
a  central  Earth  angle)  that  connects  the  generator  point  with  the  subsatellite  point,  specifically, 
from  the  magnitude  of  6e  and  its  orientation  relative  to  the  subsatellite  point  A  free  parameter 
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Application  of  the  technique  first  requires  calculation  of  two  auxiliary  angles  that  are 
derived  from  the  beam  centroid  and  satellite  positions.  These  are  the  elevation  and  azimuth 
angles,  denoted  eand  az  respectively,  of  the  beam  centroid  with  respect  to  the  satellite.  The 
angle  c  is  depicted  in  figure  3- 12  and  evaluated  from  the  dot  product  as 


e  =  arccosl - 

V  rd 

The  angle  az,  illustrated  by  the  angle  x  in  figure  3-14,  is  defined  in  the  usual  manner  as  the 
positive  angle  measured  clockwise  from  north  to  the  pointing  location.  In  general,  az  is 
calculated  as  a  function  of  the  smallest  angle  X  measured  from  north  to  the  pointing  location  as 
follows: 


f  Z,  if  ysat  lies  to  the  left  of  ycen 
[360*-£,  if  Ysat  I'es  to  the  right  of  ycen 


where  and  ycen  denote  the  longitudes  of  the  subsatellite  and  beam  centroid  points, 
respectively.  The  adjustment  that  occurs  when  ysat  lies  to  the  right  of  ycen  is  necessary  since 
in  such  cases,  az,  by  definition,  assumes  values  larger  than  180*.  The  angle  x>  in  turn,  is  a 
function  of  the  angle  b0  (also  shown  in  figure  3-14).  Both  x  and  bQ  are  calculated  from  the 
spherical  triangle  law  of  cosines  for  sides  as  follows: 


X  =  arccos 


cos(90*-SceB )  -  cos(90*-g,at  )cos  (b0 ) 
sin(90'-6sat)sin(bo) 


(3-11) 


b0  =  arccos(cos(90*-5sa/  )cos(90*-$ceB )  +  sin(90*-5sat  )sin(90*-5ceB)cos(d/ongc<!B)) 

with  Ssat  and  Scen  denoting  the  latitudes  of  the  subsatellite  and  centroid  locations,  and  Alongccn 
taken  as  the  magnitude  of  the  difference  between  Ysat  and  ycen- 

The  individual  steps  of  the  technique  are  outlined  below.  The  calculations  rely  heavily  on 
trigonometric  relationships  for  spherical  triangles.  The  details  are  deferred  to  reference  1: 

1)  Calculate /?=  (beam  width)/2 

2)  Set  the  free  parameter  y  0*  £  y<  360* 

3)  Calculate  the  side  9  and  then  the  interior  angle  0  of  the  spherical  triangle  shown  in  figure  3-13: 

6  -  arccos(cos  ecos/3  +  sin  £  sin/?  cos  y) 

0  =  arcsin(sin/Jsiny/sin0)  (3-12) 
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north  pole 


Figure  3-14.  Azimuth  Angle 

4)  Calculate  the  angle  on  the  surface  of  the  Barth,  9e,  which  corresponds  to  9 ;  assume  a 
constant  Earth  radius,  Re: 


9e  =  arcs  in 


(  rsinfl) 
1  «e  ) 


0 


(3-13) 


5)  Calculate  the  generator  point  latitude,  Sgen,  and  the  difference.  Along gen,  in  the  generator 
point  longitude  relative  to  the  subsatellite  longitude: 


Sgen  =  arcsin(sin  Ssat  cos  9e  +  cos  Ssat  sin  9e  cos  (az  +  0)) 


(3-14) 


Along gen  -  arcsin 


sin(az  +  0)sinflg 
^  cos  (Sgen)  J 


(3-15) 


6)  Increment  y  and  repeat  steps  (3)  through  (5). 

7)  Repeat  step  (6)  until  y  sweeps  out  360'. 

ESP  uses  a  fixed  step  size  for  y  of  2°. 8 125,  so  chosen  to  yield  128  generator  points  for 
every  projection.  This  value  was  found  to  produce  projections  that  are  smooth  in  appearance 
regardless  of  the  size  of  the  projection  contour. 
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In  the  above  formulation,  there  are  problems  under  certain  circumstances  with  inverse 
trigonometric  functions  returning  the  incorrect  quadrant  One  example  is  evaluation  of  the 
angle  <f>  by  equation  3-12  when  £<  f3,  i.e.,  when  the  projection  contains  the  subsatellite  point. 
The  problem  arises  because  in  this  circumstance,  0  assumes  values  larger  in  magnitude  than 
90*,  however,  the  arcsine  function  by  itself  is  incapable  of  distinguishing  angles  in  quadrants 
II  and  m  from  quadrants  I  and  IV  (by  default  it  returns  values  in  quadrant  I  when  positive  and 
quadrant  IV  when  negative).  A  similar  problem  occurs  with  equation  3-15  when  the  satellite  is 
near  the  poles.  ESP  handles  these  circumstances  by  supplementing  the  sine  evaluations  with 
cosine  evaluations,  and  using  a  “smart”  arctangent  function,  which  is  described  in  section  4.6. 
(A  similar  problem  is  generally  not  encountered  with  the  arccosine  function  since  in  many 
cases  the  angle  involved  is  known  to  fall  within  the  range  of  arccosine  values,  i.e.,  within  the 
interval  [0*,  180*].) 

The  above  technique  has  also  been  modified  so  that  it  is  consistent  with  the  Earth  model 
described  by  equation  3-6,  in  which  the  radius  of  the  Earth  is  a  function  of  geocentric  latitude. 
An  iteration  technique  is  used  on  steps  (4)  and  (5)  as  follows.  An  initial  value  8eo  is  calculated 
from  equation  3-13  with  Re  set  equal  to  the  radius  of  the  Earth  at  the  centroid  location. 

Equations  3-14  and  3-15  are  evaluated  with  6e  set  equal  to  6e  to  arrive  at  initial  generator 
point  coordinates.  The  radius  of  the  Earth  corresponding  to  the  calculated  generator  point 
latitude  is  then  used  in  equation  3-13  to  calculate  a  new  value  for  6e,  which  in  turn  is  used  in 
equations  3-14  and  3-15  to  calculate  new  generator  point  coordinates.  The  process  is  repeated 
until  subsequent  values  of  the  generator  point  coordinates  are  sufficiently  close. 

There  are  also  a  number  of  potential  singularities  (e.g.  zero  denominators)  in  the  above 
formulation  that  require  special  attention.  Some  are  easily  eliminated  by  considering  their 
geometric  inteipretation.  For  example,  equation  3-11  need  not  be  evaluated  when  b0  is  equal  to 
zero  since  this  implies  az  is  equal  to  zero.  Similarly,  equation  3-12  (and  all  subsequent 
equations  as  well)  is  not  necessary  when  6  is  equal  to  zero,  since  in  this  case,  the  generator 
point  coordinates  are  identical  to  the  subsatellite  point  coordinates.  Another  example  is 
evaluation  of  costp  (used  to  supplement  equation  3-12)  when  e  =  0.  In  this  case,  the  beam 
pointing  location  is  the  subsatellite  point  and  the  spherical  triangle  used  to  define  <f>  (shown  in 
figure  3-13)  does  not  exist.  Geometrically,  <p  becomes  a  free  parameter  and  is  taken  equal  to  y, 
thereby  eliminating  the  use  of  the  cosine  function.  In  a  few  other  instances,  namely  equation 
3-11  when  Ssat  is  equal  to  ±90*  and  equation  3-15  when  8oen  is  equal  to  ±90*,  the  singularities 
are  avoided  by  substituting  a  latitude  value  which  is  very  close,  as  described  in  section  4.6. 

This  is  done  to  preserve  the  geometry  upon  which  the  algorithm  is  based  and  thus  ensure  a 
continuous  projection  curve;  it  is  necessary  since  the  longitude  of  a  point  at  the  poles  is 
arbitrary. 

Throughout  the  algorithm,  ESP  checks  to  see  if  the  beam  shines  over  the  horizon.  This 
occurs  when  there  exist  values  of  6,  calculated  by  step  (3),  that  exceed  the  maximum  possible  6, 
denoted  ©mo*  and  shown  in  figure  3-15.  The  check  is  important  since  in  general  6  greater  than 
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Figure  3-15.  Maximum  0  Value 


Omax  means  the  planar  triangle  which  is  the  basis  for  calculating  0e  does  not  exist.  The  result  is 
that  evaluation  of  equation  3-13  requires  evaluation  of  the  arcsine  function  with  an  argument 
greater  than  one.  Since  Omax  is  a  function  of  the  radius  of  the  Earth,  which  is  variable,  an  initial 
guess  is  conservatively  chosen  as 


*  .  ( 6356.75^ 

Omax  *  arcsml — - — I 

This  approximation  is  considered  conservative  since  it  represents  the  minimum  possible  (Wx, 
corresponding  to  the  Earth’s  radius  at  the  poles.  Thus,  when  0  exceeds  Omax,  there  may  or  may 
not  be  a  problem.  In  such  cases,  two  potential  generator  points  are  calculated.  The  first  is 
calculated  in  the  usual  manner  as  described  above  with  0e  set  equal  to  an  estimate  of  the 
maximum  possible  value,  if  necessary,  to  avoid  die  arcsine  problem.  The  second  point  is 
calculated  to  fall  on  the  horizon  by  iterating  on  0e  until  the  elevation  angle  at  the  resultant 
generator  point  is  sufficiently  close  to  zero.  Subsequent  values  for  0e  are  determined  from 
previous  values  as  follows 


=  0e.  +  A 


where  the  parameter  A  is  chosen  relative  to  the  elevation  angle  as 

A  =  (elevation  angle)  /  2 

This  choice  is  based  on  an  analysis  of  the  sensitivity  of  elevation  angle  to  0e,  which  showed 
that  as  the  elevation  angle  varied  from  +0.5  to  -0.5  degrees,  0e  increased  by  an  approximate 
equal  amount  of  one  degree,  regardless  of  orbital  altitude.  In  all  of  the  test  cases,  this 
technique  worked  well,  with  the  number  of  iterations  eight  or  fewer  and  most  often  four. 

Given  the  horizon  point,  the  estimate  for  Omax  is  updated  based  on  the  corresponding  radius  of 
the  Earth.  A  choice  is  then  made  between  the  two  potential  points  according  to  whether  or  not 
6  exceeds  the  updated  value.  ESP  plots  and  connects  all  of  the  projection  points  up  to  and 
including  any  end  points  that  fall  directly  on  the  horizon.  Closure  of  a  projection  that  shines 
over  the  horizon  is  achieved  by  overlaying  a  plot  of  the  horizon  contour,  whose  generation  is 
described  in  section  3.1 1. 
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3.10  LOSS  CONTOURS 


In  the  case  in  which  a  single  beam  projection  is  of  interest  and  the  beam  width 
corresponding  to  3  dB  loss  (specified  by  the  user)  is  less  than  30*,  ESP  offers  the  option  to 
plot  contours  of  equal  loss.  A  circular  aperture  uniform  illumination  beam  pattern  function 
(one  of  the  more  common  pattern  functions  for  circular  conical  beams)  is  used.  The  following 
expression  is  used  to  approximate  voltage,  as  a  function  of  the  angle  rj  off  boresight: 


(3-16) 


where  T]0  represents  the  beam  width  corresponding  to  3  dB  loss,  J\(x)  is  the  Erst  order  Bessel 
function  of  the  first  kind,  and  K\  is  a  constant  equal  to  3.2328.  This  approximation,  which 
removes  the  explicit  dependency  of  voltage  on  the  radius  of  aperture  and  wavelength,  is 
reasonably  accurate  for  values  of  T]0  less  than  or  equal  to  30°  (which  correspond  to  a  radius  of 
aperture  less  than  or  approximately  equal  to  the  wavelength).  For  such  values,  the  3  dB  beam 
width  is  accurate  to  within  0.3  degrees,  and  the  corresponding  loss  value  is  accurate  to  within 
0.06  dB. 

ESP  takes  T]0,  by  default,  as  the  value  entered  by  the  user  in  the  beam  width  field  of  the 
dialog  box.  The  user  has  the  option  of  plotting  contours  with  either  a  uniform  degree  or  dB 
step.  In  the  case  of  the  former,  a  loop  is  performed  on  r?  with  loss  calculated  as  a  function  of 
voltage  as  given  below: 

Loss  =  201og(P(7j)) 

In  the  case  of  the  latter,  the  problem  is  solved  backward  with  a  loop  over  loss  values.  The 
result  is  that  for  a  given  loss  value,  equation  3-16  must  be  solved  for  tj.  Rearranged,  the 
problem  is  to  solve 

/(*?)  =  Jl(C2  sin  *?)  “  clc2  sin  rj  =  0  (3_17j 

where 
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Since  equation  3-17  is  transcendental  in  rj  (i.e.,  it  cannot  be  solved  explicitly  for  rj),  it  is 
solved  using  a  Newton  iteration  scheme  similar  to  the  one  described  in  section  3.5,  with  the 
function 77)  given  as  above  and  f'(tj)  given  by 


/'(*»)- 


70(C2sin  77)- 


Jl(C2sinri)y 
(C2  sin  tj)  y 


C2  cos  77  -  C]C2  cos  77 


where  7o(jc)  is  the  zero  order  Bessel  function  of  the  first  kind.  For  a  given  loss  value,  an  initial 
guess  for  77  is  taken  as  one-half  the  beam  width  of  the  previous  contour.  This  technique  was 
found  to  work  quite  well,  with  convergence  occurring  in  most  cases  within  four  or  five 
iterations. 

Regardless  of  the  step  units,  contours  are  plotted  for  loss  values  up  to  but  not  including 
-17.6  dB.  This  value  corresponds  to  the  power  at  the  first  side  lobe  of  the  beam  pattern 
function,  where  the  notion  of  a  contour  breaks  down.  The  Bessel  functions  are  evaluated 
according  to  polynomial  approximations  provided  in  reference  22. 


3.11  HORIZON  GENERATION 

A  horizon  point  is  defined  as  a  ground  location  from  which  the  elevation  angle  to  the 
satellite  is  zero.  A  horizon  contour  is  a  curve  representing  the  set  of  points  that  lie  on  the 
horizon  for  a  particular  satellite  at  a  particular  instant.  This  contour  is  often  of  interest  since  it 
defines  a  boundary  beyond  which  ground  sites  are  incapable  of  seeing  the  satellite  at  some 
point  in  time  (although  points  on  and  along  the  inner  edge  of  this  contour  are  often  unable  to 
see  the  satellite  due  to  terrain  irregularities  and  tall  buildings).  The  technique  for  generating  the 
horizon  contour  is  similar  to  the  technique  for  generating  satellite  ground  tracks  and  spot  beam 
contours  in  that  a  number  of  discrete  points  along  the  contour  are  found  and  connected.  In  the 
notation  of  section  3.8,  the  problem  is  to  find  all  right  ascension-declination  pairs,  (do,  S0), 
such  that  d  •  a  =  0  (refer  to  figure  3-11).  Substitution  of  the  vectors  d  and  a,  with  rsue 
denoting  the  vector  to  the  horizon  point  with  coordinates  (aQ  S0)  and  otherwise  the  same 
notation  as  used  in  section  3.8,  results  in  the  following  equation: 

cosao+psinao  +  Y  =  0  f3,8x 


where 


-  _  rsite{rsite  ~  bcos80  -  zsin  80) 
x(b-rsilecos80) 
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Fora  specific  value  for  So,  the  above  coefficients,  ft  and  y,  are  constant.  Elimination  of 
sinab  in  equation  3-18  reduces  the  problem  to  solving  a  quadratic  in  cosab: 


acos2(a0 ) + b  cos  a0  +  c  =  0 


(3-19) 


where 


a  =  l  +  p2 


b  =  2y 


c-y2  -p2 

The  idea  is  then  to  set  80  and  use  the  quadratic  formula  to  solve  equation  3-19  for  cosab- 
There  may  be  at  most  two  solutions  to  the  quadratic  which  translates  into  four  possible  right 
ascension  values,  since  a  solution  A  for  cosab  implies  that  both  arccosA  and  360*  -  arccosA  are 
possibilities.  The  elevation  angle  for  each  possible  coordinate  pair  is  then  calculated  as 
described  in  section  3.8  to  determine  which  possibilities  correspond  to  points  along  the 
horizon.  A  complete  set  of  horizon  points  is  obtained  by  iterating  this  procedure  over  latitude 
values  S0  ranging  from  -90*  to  90*.  Right  ascension  coordinates  are  convened  to  longitude 
values  for  plotting  purposes  as  described  in  section  3.6. 
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SECTION  4 

PROGRAMMING  CONSIDERATIONS 


4.1  INTRODUCTION 

This  section  addresses  testing  and  questions  of  numerical  accuracy  within  ESP.  As  with 
any  computation-intensive  program,  measures  had  to  be  taken  to  deal  with  the  unavoidable 
errors  resulting  from  the  finite  precision  with  which  computers  represent  floating-point 
numbers. 


4.2  FLOATING  POINT  PRECISION 

ESP  uses  double  precision  floating  point  numbers  for  time,  date,  prime  meridian,  and 
certain  horizon  calculations,  and  single  precision  for  all  other  floating  point  operations.  ESP 
uses  Julian  dates,  where  the  basic  unit  of  time  is  a  day.  Julian  date  values  are  currently  so 
large  (for  example,  1  January  1990  is  2447892.5)  that  single  precision  floating  point  numbers 
give  a  resolution  of  about  one  day.  Double  precision  is  used  to  represent  hours,  minutes,  and 
seconds.  This  is  particularly  important  when  calculating  the  prime  meridian  location  and  the 
orbital  elements  of  the  Sun  as  functions  of  time.  In  both  examples,  double  precision  is  used  to 
ensure  sufficient  precision  in  determining  the  time  elapsed  from  the  1950  reference  epoch  to  the 
present. 


4.3  CONVERGENCE  AND  NUMERICAL  TOLERANCE 

ESP  performs  several  operations  that  require  the  difference  between  two  values  to  be 
“close.”  An  example  is  the  use  of  Newton’s  method,  which  derives  successive  estimates  of 
the  root  of  an  equation  from  previous  estimates.  ESP  uses  Newton’s  method  to  solve  Kepler’s 
equation  for  eccentric  anomaly  and  to  convert  a  spot  beam  decibal  loss  value  into  beam  width. 
For  termination,  this  algorithm  requires  a  maximum  tolerance  value  for  the  difference  between 
successive  estimates. 

In  general,  it  is  possible  for  the  difference  between  estimates  never  to  fall  below  a  given 
tolerance.  If  this  happens  in  an  iterative  algorithm,  the  program  will  fall  into  an  infinite  loop. 
At  the  same  time,  it  is  important  that  the  tolerance  be  small  enough  to  insure  sufficient  precision 
in  the  result.  It  is  therefore  important  to  allow  a  tolerance  slightly  larger  than  the  minimum 
non-zero  difference  between  any  two  estimates.  (Unlike  in  real  arithmetic,  in  floating  point 
arithmetic  there  can  be  a  minimum  non-zero  difference  between  values.)  This  minimum  will 
vary  with  the  size  of  the  estimates.  The  program  therefore  chooses  tolerances  relative  to  the 
estimates,  using  the  formula 

(l£ll  +  IE2l)(*) 
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where  E\  and  £2  are  the  two  values  being  compared,  and  k  is  10*7  for  single  precision  and 
2x1 for  double  precision.  This  formula  results  in  tolerances  very  near  the  limits  of  the 
Macintosh’s  precision. 

In  stand-alone  tests  of  the  Kepler  function,  estimates  invariably  converged  to  a  difference 
of  zero  within  five  iterations.  However,  since  there  has  not  been  a  rigorous  proof  that  the 
difference  between  successive  estimates  will  always  converge  to  zero,  the  tolerance  is  an 
important  safeguard  to  prevent  an  infinite  loop.  In  similar  tests  of  the  spot  beam  decibel  loss- 
to-angle  conversion,  the  tolerance  was  found  to  be  necessary,  since  the  difference  between 
estimates  did  not  always  converge  to  exactly  zero  (although  the  difference  invariably  converged 
to  within  the  tolerance  within  seven  iterations,  and  usually  within  five). 

In  the  unlikely  event  that  convergence  does  not  take  place  after  a  reasonable  number  of 
iterations,  an  error  message  is  presented  to  the  user,  and  calculation  proceeds  normally  after  the 
user  dismisses  the  error  message.  This  error  message  has  never  been  generated  in  actual  use 
within  D063.  “Reasonable”  is  defined  as  ten  iterations  in  this  case. 

A  tolerance  is  also  used  to  calculate  the  horizon.  This  calculation  uses  the  quadratic 
formula  to  locate  points  along  the  horizon.  In  particular,  the  discriminant  of  the  formula 
determines  whether  or  not  such  points  exist  at  a  given  latitude.  If  the  absolute  value  of  the 
discriminant  falls  below  the  tolerance  calculated  from  the  formula  above  (with  £1  equal  to  zero 
and  £2  equal  to  the  discriminant),  the  discriminant  is  assumed  to  be  zero.  This  allows  certain 
points  to  be  located  which  would  otherwise  be  skipped  due  to  a  very  small  negative 
discriminant.  Since  the  program  searches  for  points  at  every  degree  of  latitude,  skipping  such 
a  point  would  create  an  error  in  the  graphical  output 


4.4  THE  RUNGE-KUTTA  PROCEDURE 

ESP  uses  two  levels  of  time  steps  in  its  calculations.  The  higher  level  time  step  is  the  one 
entered  by  the  user  in  the  dialog  box.  In  order  to  ensure  accurate  results  from  the  numerical 
integration,  the  calculations  are  done  at  intervals  of  about  one  minute  or  less.  These  one- 
minute  intervals  will  be  referred  to  as  the  internal  time  interval;  the  one-minute  step  size  will  be 
referred  to  as  the  internal  time  step.  It  is  at  the  internal  time  intervals  that  ESP  calculates  the 
satellite  position  and  velocity.  (The  internal  time  step  is  derived  by  dividing  the  user-specified 
time  step  into  equal  intervals.  Therefore,  the  internal  time  step  wifi  not  be  exactly  one  minute  if 
the  user-specified  step  size  is  not  an  exact  integer  number  of  minutes.) 

Each  call  to  the  Runge-Kutta  procedure  calculates  the  satellite  position  and  velocity  over 
one  user-specified  (higher  level)  time  interval.  The  procedure  itself  calculates  the  satellite 
position  and  velocity  many  times,  using  the  internal  time  step.  If  each  new  time  were  derived 
by  adding  the  internal  time  step  to  the  previous  time,  there  would  be  the  potential  for  significant 
round-off  errors.  The  time  at  die  end  of  the  procedure’s  calculations  might  not  match  the  time 
that  would  be  expected  based  on  the  user  inputs.  More  seriously,  the  output  would  show 
satellite  positions  and  velocities  that  do  not  match  the  printed  times. 
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Although  the  Runge-Kutta  procedure  can  take  measures  to  minimize  any  mismatch  between 
the  calculated  time  at  the  end  of  its  run  and  the  actual  user  specified  time,  it  cannot  guarantee 
that  the  two  numbers  will  be  equal  each  time.  Therefore,  each  time  the  procedure  is  called 
(except  for  the  first),  the  calculated  ending  time  from  the  previous  call  is  used  as  the  starting 
time  for  the  new  higher  level  interval.  This  is  done  in  order  to  insure  consistency  between  the 
time  values  and  the  satellite  positions  and  velocities  from  one  call  of  the  procedure  to  the  next. 
The  user-specified  interval  ending  time  is  passed  to  the  procedure,  giving  it  a  target  time  and 
insuring  that  the  times  used  in  the  calculations  do  not  stray  from  the  user-specified  times. 

Within  the  Runge-Kutta  procedure  itself,  measures  were  taken  to  minimize  round-off 
errors  in  the  internal  time  interval  calculation.  The  procedure  determines  the  internal  time 
values  in  a  two  step  process,  relying  entirely  on  integers  for  loop  control.  First,  the  number 
of  iterations  requirwi  over  the  higher  level  interval  is  determined.  Then,  on  every  pass  through 
the  loop,  the  precise  values  of  the  beginning  and  the  end  of  each  internal  time  interval  are 
recalculated  from  the  loop  control  values  and  the  initial  time  of  the  higher  level  interval  passed 
into  the  procedure.  The  variable  that  contains  the  initial  time  of  the  higher  level  interval  is 
never  changed  during  this  process.  In  this  way,  the  cumulative  errors  of  repeated  floating 
point  additions  are  avoided.  Upon  exit  from  the  procedure,  the  ending  time  of  the  higher  level 
interval  is  calculated  in  the  same  manner  based  on  its  initial  time.  Only  one  floating  point 
addition  is  used  for  the  ending  time  value. 

The  Runge-Kutta  algorithm  inherently  introduces  error,  since  the  equations  are  derived 
from  truncated  Taylor  series  expansions  where  terms  of  the  order  (Ar)5  and  higher  are 
disregarded.  As  a  means  of  assessing  the  magnitude  of  the  error  due  to  truncation  and  a  fixed 
time  step  At,  results  for  a  variety  of  orbit  types  (low-Earth,  Molniya,  and  geosynchronous) 
were  compared  with  results  from  ASAP,  which  is  significantly  more  sophisticated  in  its 
implementation  of  an  eighth-order  algorithm  with  variable  step  size  control.  The  outcome  was 
quite  favorable;  over  a  24-hour  time  interval,  position  vector  components  agreed  to  within  one 
kilometer  at  worst  and  often  within  tens  of  meters.  This  consistency  of  results  is  a  source  of 
confidence  as  well  as  validation,  since  ASAP  is  in  turn  a  well-tested  tool  with  results  in 
agreement  with  other  orb’ t  propagators  and  more  importantly,  actual  tracking  data. 

4.5  THE  MAIN  LOOPS 

For  the  main  loops  over  time,  longitude  and  latitude,  all  of  the  loop  control  parameters  are 
floating  point.  Due  to  the  inexact  nature  of  floating  point  control  parameters,  it  is  not  unusual 
for  a  loop  to  iterate  one  time  too  few  when  using  such  parameters.  This  can  be  demonstrated 
with  an  example  from  a  hypothetical  machine  with  a  precision  of  three  decimal  digits.  The 
number  two-thirds  would  be  represented  as  0.667.  A  loop  from  0.0  to  10.0  by  two-thirds 
would  execute  only  fifteen  times  if  10.0  /  0.667  is  truncated,  and  would  end  at  9.33  instead  of 
10.0. 

To  solve  this  problem,  the  correct  number  of  iterations  is  determined  before  any  loop  is 
executed.  A  tolerance  factor  relative  to  the  size  of  the  difference  between  the  stopping  and 
starting  parameters  is  added  to  that  difference.  The  adjusted  difference  is  then  divided  by  the 
step  size.  When  the  loop  is  executed,  an  integer  variable  is  used  for  counting,  using  the 
derived  number  of  iterations.  The  proper  floating  point  value  for  each  pass  is  derived  in  turn. 
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The  tolerance  used  is  large  enough  to  account  for  round-off  error,  but  small  enough  to  avoid 
overstepping  the  upper  limit,  so  that,  for  example,  a  loop  from  0  to  10  with  a  step  of  1 .5 
would  not  end  at  10.5. 

If  the  range  of  the  loop  is  not  a  multiple  of  the  step  size  (within  the  tolerance),  the 
computed  number  of  iterations  is  increased  by  one  in  order  to  include  the  upper  bound  entered 
by  the  user. 


4.6  NUMERIC  ERROR  AVOIDANCE 

ESP  uses  a  “smart”  arctangent  function,  which  has  a  range  of  (-K,  it]  and  handles  90’  and 
270*  properly.  The  function  takes  two  arguments  and  is  defined  as 

arctan2  (y,  x)  =  arctan  (y/x) 

with  special  cases  to  handle  x  =  0  and  testing  of  the  signs  of  x  and  y  to  determine  the  proper 
quadrant  of  the  result  For  the  reader  familiar  with  the  FORTRAN  programming  language, 
this  function  is  equivalent  to  FORTRAN’S  function  ATAN2. 

The  program’s  arcsine  and  arccosine  implementations  use  arctan2  rather  than  the  standard 
Pascal  arctan,  so  arcsin  (1)  and  arccos  (0)  do  not  result  in  runtime  errors.  Occasionally,  as  a 
result  of  round-off  errors,  a  number  slightly  greater  than  1  or  less  than  -1  may  be  passed  to  one 
of  these  two  routines  during  the  course  of  calculation.  In  these  rare  cases  a  message  is 
presented  to  the  user,  a  1  (or  -1)  is  substituted  for  the  incorrect  argument,  and  calculation 
proceeds  normally. 

As  noted  in  section  3.9,  a  singularity  occurs  during  spot  beam  calculations  when  the 
subsatellite  point  or  a  spot  beam  generator  point  is  at  a  pole.  In  addition,  the  geometry  breaks 
down  when  either  point  is  very  near  a  pole.  In  these  instances,  89.999*  is  substituted  for  90°. 
The  substitution  does  not  affect  the  output  because  the  output  is  given  to  two  decimal  places. 


4.7  FUNCTION  AND  PROCEDURE  TESTING 

Where  possible,  the  functions  and  procedures  in  ESP  were  individually  tested  and 
validated.  Some  specific  examples  are  cited  below: 

•  The  Julian  date  function  was  verified  by  comparison  with  the  JPL  software  tool 
QUICK  [20].  This  interactive  calculator-like  program  has  been  extensively  used  and 
tested  at  JPL  since  its  initial  version  in  1971. 

•  As  described  above,  the  Newton  iteration  schemes  used  to  solve  Kepler’s  equation 
and  to  calculate  beam  width  as  a  function  of  decibel  loss  were  tested  for  convergence. 
In  addition,  the  solutions  were  verified  for  a  variety  of  examples  by  substituting  them 
back  into  the  original  equations. 
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•  The  Runge-Kutta  algorithm  was  validated  via  trigonometric  examples.  Values  for  the 
sine  and  cosine  functions  were  estimated  and  found  to  agree  with  die  actual  values  to  a 
minimum  of  five  (and  often  six  or  seven)  decimal  places. 

•  The  inverse  trigonometric  functions  were  tested  with  angles  ranging  from  0*  to  360*  to 
verify  that  the  results  returned  were  in  the  proper  quadrant 

•  Various  altitudes  were  passed  into  the  atmospheric  density  function  to  insure  it 
produced  the  same  values  as  in  the  density  input  tables. 

•  Sanity  checks  were  made  on  the  elevation  and  outage  outputs  to  insure  that  they 
reflected  the  calculated  satellite  positions. 

•  The  coded  polynomial  approximations  for  Jq(x)  and  J\  (x)  were  verified  by 
comparison  with  tabulated  values  provided  in  reference  22. 

•  The  spot  beam  projection  and  horizon  algorithms  were  tested  by  examining  numerous 
possible  scenarios  (e.g.,  for  the  relative  locations  of  the  pointing  location  and 
subsatellite  point  and  a  variety  of  orbit  types).  Doing  so  highlighted  a  number  of 
singularities  (particularly  in  the  case  of  polar  orbits)  that  have  since  been  eliminated. 
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APPENDIX 


SAMPLE  RUNS 


A.l  SUMMARY 

This  appendix  contains  a  sample  of  the  program's  outputs.  It  is  organized  into  three 
subsections  according  to  orbit  type.  A  summary  of  the  content  of  the  individual  subsections 
including  a  brief  description  of  the  orbits  presented  is  provided  below.  Unless  otherwise 
stated,  satellite  orbits  were  propagated  over  24  hours,  and  outage  feature  runs  have  the 
minimum  elevation  angle  parameter  set  to  10*  and  the  minimum  number  of  satellites  in  view 
parameter  set  to  one.  Additional  examples  which  show  how  the  program  has  been  used  since 
its  initial  release  in  1991  to  support  a  variety  of  MITRE  projects  can  be  found  in  reference  7. 

Section  A.2  provides  all  types  of  program  outputs  for  geosynchronous  orbits;  it  consists  of 
the  echoed  input  data  as  well  as  both  tabular  and  graphical  results.  It  should  be  noted  that  the 
tabular  output  shown  is  just  a  one  page  sample  in  all  cases,  since  the  complete  output  is  often 
quite  lengthy  and  would  require  numerous  pages  to  duplicate.  By  definition,  a 
geosynchronous  orbit  is  a  circular  orbit  with  period  equal  to  the  period  of  the  Earth's  rotation 
about  its  axis  (approximately  23.93  hours).  The  semimajor  axis  is  calculated  to  be  42,163  km 
from  equation  3-3.  A  satellite  in  geosynchronous  orbit  is  often  considered  ideal  for 
communications  purposes  since  it  resides  over  the  same  general  region  of  the  Earth.  This  is 
obvious  finom  the  ground  track  and  becomes  remarkably  more  pronounced  as  the  inclination 
angle  approaches  zero.  In  all  of  the  examples  in  this  section,  satellite  orbits  have  a  65* 
inclination.  Figure  A-l  shows  the  ground  track  for  a  satellite  centered  on  99*  W  longitude. 
Table  A-4  presents  hourly  elevation  angle,  range,  and  range  rate  calculations  for  the  same 
satellite  as  it  travels  through  its  orbit  and  a  set  of  four  ground  sites  located  in  western  Europe. 
Figure  A-2  is  an  outage  plot  for  a  constellation  of  three  satellites  with  parameters  Q,  a,  and  z 
chosen  so  that  each  follows  approximately  the  same  ground  track  depicted  in  figure  A-l.  Such 
a  constellation  provides  near  perfect  coverage  for  one  side  of  the  globe  (as  can  be  seen  by  the 
areas  that  have  no  markings  in  the  plot)  and  no  coverage  for  the  other  side  (as  can  be  seen  by 
the  large  region  that  falls  ir  t.  -ne  6).  The  output  time  step  for  both  the  ground  track  and 
outage  zones  is  chosen  rehj •  >'c  o  the  orbit  period  and  is  one-half  hour,  this  ensures  both  a 
smooth  ground  track  and  red  uably  accurate  coverage  results.  The  last  two  examples  show 
the  program's  spot  beam  and  loss  contours  features  for  the  same  satellite  as  in  figure  A-l, 
although  the  projections  are  calculated  at  different  points  along  its  orbit.  In  figure  A-3,  a  5’ 
beam  is  pointing  in  the  Boston  area,  whereas  in  figure  A-4,  a  4°  beam  is  pointing  directly 
below  the  satellite.  Both  plots  have  been  enlarged  via  the  zoom  feature  described  in  section 
2.4.5  under  Draw  a  Map. 

The  remaining  sections  contain  only  the  graphical  outputs.  Section  A.3  is  based  on  typical 
MOBILSATCOM  orbits  [19],  which  are  highly  inclined,  circular  low-Earth  orbits.  The 
satellite  altitude  is  approximately  1000  km  (which  corresponds  to  a  semimajor  axis  of  7378 
km)  and  the  orbital  inclination  is  80*.  Figure  A-5  shows  the  ground  track.  As  can  be  seen 
from  the  plot,  the  satellite  travels  around  the  Earth  several  times  over  the  course  of  a  day  (13.7 
to  be  exact),  since  the  orbit  period  is  small  (105  minutes)  relative  to  the  period  of  the  Earth's 
daily  rotation.  Figure  A-6  shows  the  coverage  plot  for  a  96-satellite  constellation  that  consists 
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of  eight  different  rings  equally  spaced  in  12  with  12  equally  spaced  satellites  per  ring.  Satellites 
in  adjacent  planes  are  half  way  out  of  phase,  and  a  minimum  of  two  satellites  must  be  in  view 
for  ground  sites  to  be  considered  covered.  In  both  examples  a  time  step  of  three  minutes  was 
used  in  order  to  obtain  sufficient  resolution  for  the  ground  track  and  outage  zones. 

Examination  of  the  program's  text  output  indicates  that  the  maximum  cumulative  outage  over 
the  course  of  a  day  is  4.3  hours.  An  example  of  multiple  beam  projections  is  shown  in  figure 
A-7,  in  which  there  are  three  35*  spot  beams  located  on  the  same  satellite  as  in  figure  A-5. 
Since  the  beams  shine  partially  over  the  horizon,  the  horizon  contour  is  plotted  as  well. 

Section  A.4  presents  results  for  Molniya  orbits.  These  are  defined  as  highly  eccentric 
12-hour  orbits  (semimajor  axis  equal  to  26,610  km)  with  inclination  given  by 

J_ 

V5 

This  value  is  referred  to  as  the  critical  inclination.  At  this  inclination,  the  line  of  apsides  will  on 
the  average  remain  constant  over  time.  In  contrast,  the  line  of  apsides  will  regress  (advance) 
when  the  orbit  inclination  is  greater  (less)  than  this  value  due  to  Jj  effects.  Figure  A-8  shows 
the  ground  track  for  a  satellite  in  a  Molniya  orbit  with  an  eccentricity  of  0.72  and  with  perigee 
in  the  southern  hemisphere.  In  order  to  close  the  ground  track,  it  was  necessary  to  propagate 
the  orbit  over  24  hours,  8  minutes,  and  57  seconds,  since  it  takes  half  that  time  to  generate  the 
portion  of  the  track  from  the  start  of  the  loop  centered  on  100*  W  longitude  to  the  start  of  the 
loop  centered  on  78*  E  longitude.  A  major  advantage  of  such  orbits  is  that  a  large  portion  of 
the  Earth’s  surface  is  visible  at  high  elevation  angles  for  long  intervals  of  time.  This  is  a  result 
of  the  fact  that  the  satellite  spends  a  large  percentage  of  its  time  near  apogee  (approximately  10 
hours  to  form  each  of  the  ground  track  loops  in  figure  A-8)  and  thus  a  relatively  small  period 
of  time  near  perigee.  Figure  A-9  shows  the  coverage  plot  for  a  constellation  of  three  satellites 
and  a  minimum  elevation  angle  of  30*.  All  satellites  have  perigee  in  the  southern  hemisphere 
and  £2  and  Tare  chosen  so  that  all  three  follow  approximately  the  same  ground  track  shown  in 
figure  A-8.  In  both  examples,  a  small  output  time  step  was  used  (3.6  minutes)  since  the 
satellite  is  moving  so  rapidly  through  perigee  (approximately  9.6  km/s  versus  1.6  km/s  at 
apogee). 


arccos 


66 


Table  A-l.  Geosynchronous  Orbit  Ground 
Track  -  Echoed  Inputs 


Satellite  Information 
Semimajor  axis 
Eccentricity 
Inclination 

Longitude  of  the  node 
Argument  of  perigee 
lime  of  perigee 


a  =42163.0  km 
e  =0.0 
i  =  65.0°  north 
.0=90.0°  east 
to  =  0.0°  east 
T  =  Jan  1, 1991  0:00:00 


The  time  of  the  osculating  values  is  Jan  1, 1991 0:00:00. 


Satellite  mass  1500.0  kg 

Effective  area  for  drag  8.0  square  m 

Coefficient  of  drag  2.0 

Effective  area  for  solar  radiation  10.0  square  m 
Absoiption/reflectivity  constant  1.5 


Inputs  For -Ihis-Rm 

Starting  time  of  the  run 
Ending  time  of  the  run 
Time  step 


Jan  1, 1991  0:00:00 
Jan  2, 1991  0:00:00 
0.5  hour 
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Table  A-2.  Geosynchronous  Orbit  Ground  Track  .  Sample  Tabular  Output 


Date 

Subsatellite 

Subsatellite 

Satellite 

and 

longitude 

latitude 

altitude 

lime 

(degrees! 

(degrees! 

(km! 

Jan  1, 1991 

0:00:00 

-9.41 

0.00 

35784.86 

Jan  1, 1991 

0:30:00 

-13.74 

6.81 

35785.15 

Jan  1. 1991 

1:00:00 

-17.98 

13.60 

35785.99 

Jan  1, 1991 

1:30:00 

-22.02 

20.35 

35787.33 

Jan  1, 1991 

2:00:00 

-25.74 

27.02 

35789.09 

Jan  1, 1991 

2:30:00 

-28.99 

33.58 

35791.12 

Jan  1, 1991 

3:00:00 

-31.54 

39.96 

35793.29 

Jan  1, 1991 

3:30:00 

-33.08 

46.09 

35795.44 

Jan  1, 1991 

4:00:00 

-33.19 

51.83 

35797.43 

Jan  1, 1991 

4:30:00 

-31.26 

56.98 

35799.10 

Jan  1. 1991 

5:00:00 

-26.62 

61.20 

35800.36 

Jan  1, 1991 

5:30:00 

-18.94 

64.03 

35801.10 

Jan  1, 1991 

6:00:00 

-9.07 

65.00 

35801.28 

Jan  1,  1991 

6:30:00 

0.72 

63.89 

35800.88 

Jan  1, 1991 

7:00:00 

8.20 

60.95 

35799.91 

Jan  1, 1991 

7:30:00 

12.64 

56.66 

35798.43 

Jan  1, 1991 

8:00:00 

14.42 

51.46 

35796.54 

Jan  1,  1991 

8:30:00 

14.20 

45.69 

35794.38 

Jan  1, 1991 

9:00:00 

12.57 

39.54 

35792.07 

Jan  1, 1991 

9:30:00 

9.97 

33.14 

35789.78 

Jan  1, 1991 

10:00:00 

6.68 

26.57 

35787.65 

Jan  1,  1991 

10:30:00 

2.94 

19.90 

35785.84 

Jan  1,  1991 

11:00:00 

-1.12 

13.15 

35784.45 

Jan  1,  1991 

11:30:00 

-5.37 

6.35 

35783.59 

Jan  1,  1991 

12:00:00 

-9.70 

-0.46 

35783.31 

Jan  1, 1991 

12:30:00 

-14.02 

-7.27 

35783.64 

Jan  1,  1991 

13:00:00 

-18.24 

-14.06 

35784.54 

Jan  1, 1991 

13:30:00 

-22.27 

-20.80 

35785.96 

Jan  1, 1991 

14:00:00 

-25.96 

-27.47 

35787.80 

Jan  1,  1991 

14:30:00 

-29.17 

-34.01 

35789.93 

Jan  1, 1991 

15:00:00 

-31.66 

-40.38 

35792.22 

Jan  1, 1991 

15:30:00 

-33.13 

-46.49 

35794.49 

Jan  1,  1991 

16:00:00 

-33.12 

-52.20 

35796.61 

Jan  1, 1991 

16:30:00 

-31.02 

-57.30 

35798.42 

Jan  1, 1991 

17:00:00 

-26.17 

-61.44 

35799.82 

Jan  1, 1991 

17:30:00 

-18.30 

-64.16 

35800.71 

Jan  1, 1991 

18:00:00 

-8.35 

-64.99 

35801.03 

Jan  1, 1991 

18:30:00 

1.33 

-63.75 

35800.78 

Jan  1,  1991 

19:00:00 

8.61 

-60.70 

35799.95 

Jan  1, 1991 

19:30:00 

12.86 

-56.33 

35798.65 

Jan  1, 1991 

20:00:00 

14.48 

-51.09 

35796.93 

Jan  1,  1991 

20:30:00 

14.14 

-45.28 

35794.91 

Jan  1,  1991 

21:00:00 

12.44 

-39.11 

35792.76 

Jan  1, 1991 

21:30:00 

9.78 

-32.70 

35790.63 

Jan  1,  1991 

22:00:00 

6.46 

-26.13 

35788.64 

Jan  1,  1991 

22:30:00 

2.68 

-19.44 

35786.95 

Jan  1,  1991 

23:00:00 

-1.39 

-12.69 

35785.69 

Jan  1,  1991 

23:30:00 

-5.64 

-5.89 

35784.93 
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Figure  A-l.  Geosynchronous  Orbit  Ground  Track  -  Graphical  Output 


Table  A-3.  Geosynchronous  Orbit  Elevation 
Angles  -  Echoed  Inputs 


Satellite  Motmapon 


Semimajor  axis 

Eccentricity 

Inclination 

Longitude  of  the  node 
Argument  of  perigee 
Time  of  perigee 


a  =42163.0  km 
e  =  0.0 
i  =  65.0°  north 
C2  -  90.0°  east 
co  =  0.0°  east 
T  =  Jan  1, 19910:00:00 


The  time  of  the  osculating  values  is  Jan  1, 1991  0:00:00. 


Satellite  mass  1500.0  kg 

E 1  tiective  area  for  drag  8.0  square  m 

Coefficient  of  drag  2.0 

Effective  area  for  solar  radiation  1 0.0  square  m 
Absorption/reflectivity  constant  1 . 5 


Inputs  For  This  Run 

Latitude  of  first  site  45.0°  north 

Latitude  of  last  site  50.0°  north 

Latitude  step  size  5.0°  north 


Longitude  of  first  site 
Longitude  of  last  site 
Longitude  step  size 

Starting  time  of  the  run 
Ending  time  of  the  run 
Time  step 


0.0°  east 
10.0°  east 
10.0°  east 

Jan  1, 1991  0:00:00 
Jan  2, 1991  0:00:00 
1.0  hour 
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Tabic  A-4.  Geosynchronous  Orbit  Elevation  Angles  -  Sample  Tabular  Output 


Date 

and 

Site 

longitude 

Site 

latitude 

Elevation 

angle 

Range 

Range 

rate 

time 

(degrees! 

(degrees) 

(degrees) 

(km) 

(km/hr) 

Jan  1, 1991 

0:00:00 

0.00 

45.00 

37.15 

37996.08 

-1064.89 

Jan  1, 1991 

0:00:00 

0.00 

50.00 

31.81 

38442.80 

-1161.25 

Jan  1, 1991 

0:00:00 

10.00 

45.00 

34.52 

38212.26 

-931.75 

Jan  1, 1991 

0:00:00 

10.00 

50.00 

29.62 

38637.02 

-1041.15 

Jan  1, 1991 

1:00:00 

0.00 

45.00 

49.29 

37120.71 

-682.18 

Jan  1, 1991 

1:00:00 

0.00 

50.00 

44.41 

37450.37 

-817.47 

Jan  1, 1991 

1:00:00 

10.00 

45.00 

44.33 

37457.39 

-580.54 

Jan  1. 1991 

1:00:00 

10.00 

50.00 

40.31 

37753.78 

-724.79 

Jan  1, 1991 

2:00:00 

0.00 

45.00 

58.03 

36620.92 

-331.09 

Jan  1, 1991 

2:00:00 

0.00 

50.00 

54.63 

36804.07 

-484.46 

Jan  1, 1991 

2:00:00 

10.00 

45.00 

50.83 

37030.08 

-291.42 

Jan  1,  1991 

2:00:00 

10.00 

50.00 

48.56 

37174.26 

-447.04 

Jan  1, 1991 

3:00:00 

0.00 

45.00 

62.34 

36419.11 

-96.42 

Jan  1, 1991 

3:00:00 

0.00 

50.00 

61.62 

36451.44 

-239.92 

Jan  1, 1991 

3:00:00 

10.00 

45.00 

54.25 

36831.39 

-129.83 

Jan  1, 1991 

3:00:00 

10.00 

50.00 

54.35 

36825.97 

-268.82 

Jan  1, 1991 

4:00:00 

0.00 

45.00 

63.36 

36378.05 

-9.77 

Jan  1,  1991 

4:00:00 

0.00 

50.00 

65.66 

36283.13 

-116.84 

Jan  1,  1991 

4:00:00 

10.00 

45.00 

56.26 

36723.95 

-105.83 

Jan  1,  1991 

4:00:00 

10.00 

50.00 

58.67 

36598.43 

-203.52 

Jan  1, 1991 

5:00:00 

0.00 

45.00 

63.87 

36358.62 

-44.51 

Jan  1, 1991 

5:00:00 

0.00 

50.00 

68.40 

36183.15 

-96.20 

Jan  1, 1991 

5:00:00 

10.00 

45.00 

58.92 

36587.77 

-176.54 

Jan  1,  1991 

5:00:00 

10.00 

50.00 

63.10 

36392.46 

-216.53 

Jan  1, 1991 

6:00:00 

0.00 

45.00 

66.03 

36271.47 

-131.30 

Jan  1,  1991 

6:00:00 

0.00 

50.00 

71.71 

36076.04 

-119.03 

Jan  1, 1991 

6:00:00 

10.00 

45.00 

63.74 

36365.07 

-264.25 

Jan  1,  1991 

6:00:00 

10.00 

50.00 

69.05 

36161.57 

-240.56 

Jan  1,  1991 

7:00:00 

0.00 

45.00 

70.62 

36108.01 

-182.38 

Jan  1,  1991 

7:00:00 

0.00 

50.00 

76.22 

35956.23 

-108.72 

Jan  1, 1991 

7:00:00 

10.00 

45.00 

71.42 

36083.16 

-280.81 

Jan  1, 1991 

7:00:00 

10.00 

50.00 

77.25 

35933.55 

-198.49 

Jan  1,  1991 

8:00:00 

0.00 

45.00 

76.54 

35945.58 

-118.02 

Jan  1, 1991 

8:00:00 

0.00 

50.00 

79.19 

35892.90 

3.34 

Jan  1, 1991 

8:00:00 

10.00 

45.00 

81.82 

35851.60 

-154.36 

Jan  1, 1991 

8:00:00 

10.00 

50.00 

86.37 

35807.38 

-29.44 

Jan  1, 1991 

9:00:00 

0.00 

45.00 

77.25 

35926.00 

106.47 

Jan  1, 1991 

9:00:00 

0.00 

50.00 

73.75 

36009.08 

252.12 

Jan  1, 1991 

9:00:00 

10.00 

45.00 

83.01 

35832.25 

143.35 

Jan  1, 1991 

9:00:00 

10.00 

50.00 

77.33 

35924.09 

285.89 

Jan  1, 1991 

10:00:00 

0.00 

45.00 

67.31 

36208.38 

478.33 

Jan  1, 1991 

10:00:00 

0.00 

50.00 

61.76 

36437.41 

619.82 

Jan  1, 1991 

10:00:00 

10.00 

45.00 

67.97 

36184.36 

577.07 

Jan  1,  1991 

10:00:00 

10.00 

50.00 

62.24 

36415.72 

709.07 

Jan  1,  1991 

11:00:00 

0.00 

45.00 

52.72 

36907.14 

922.26 

Jan  1, 1991 

11:00:00 

0.00 

50.00 

47.07 

37264.76 

1034.54 

Jan  1,  1991 

11:00:00 

10.00 

45.00 

51.18 

37000.12 

1051.06 

Jan  1, 1991 

11:00:00 

10.00 

50.00 

45.84 

37348.47 

1150.27 
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Table  A-5.  Geosynchronous  Orbit  Outage 
Zones  •  Echoed  Inputs 


Satellite  Number  1 

Semimajor  axis 

Eccentricity 

Inclination 

Longitude  of  the  node 
Argument  of  perigee 
Time  of  perigee 


a  =42163.0  km 
e  =0.0 
i  =65.0°  north 
f2=  90.0°  east 
co  =  0.0°  east 
t  =  Jan  1, 1991  0:00:00 


The  time  of  die  osculating  values  is  Jan  1, 1991  0:00:00. 


Satellite  mass 
Effective  area  for  drag 
Coefficient  of  drag 
Effective  area  for  solar  radiation 
Absorption/reflectivity  constant 


1500.0  kg 
8.0  square  m 
2.0 

10.0  square  m 
1.5 


Satellite  Number  2 

Semi-major  axis 

Eccentricity 

Inclination 

Longitude  of  the  node 
Argument  of  perigee 
Time  of  perigee 


a  =42163.0  km 
e  =0.0 
i  =65.0°  north 
i2  =  210.0°  east 
co  =  240.0°  east 
x  =Jan  1,1991  0:00:00 


The  time  of  the  osculating  values  is  Jan  1, 1991  0:00:00. 


Satellite  mass 
Effective  area  for  drag 
Coefficient  of  drag 
Effective  area  for  solar  radiation 
Absoiption/reflectivity  constant 


1500.0  kg 
8.0  square  m 
2.0 

10.0  square  m 
1.5 


Satellite  Number  3 

Semi-major  axis 

Eccentricity 

Inclination 

Longitude  of  the  node 
Argument  of  perigee 
Time  of  perigee 


a  =42163.0  km 
e  =0.0 
i  =65.0°  north 
J2=  330.0°  east 
co  =  120.0°  east 
x  =  Jan  1,19910:00:00 


The  time  of  the  osculating  values  is  Jan  1, 1991  0:00:00. 


Satellite  mass 
Effective  area  for  drag 
Coefficient  of  drag 
Effective  area  for  solar  radiation 
Absorption/reflectivity  constant 


1500.0  kg 
8.0  square  m 
2.0 

10.0  square  m 
1.5 
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Table  A-5.  (Concluded) 


Inputs.  For  This  Run 


Latitude  of  first  site 

*90.0°  north 

Latitude  of  last  site 

90.0°  north 

Latitude  step  size 

5.0°  north 

Longitude  of  first  site 

-180.0°  east 

Longitude  of  last  site 

180.0°  east 

Longitude  step  size 

5.0°  east 

Starting  time  of  the  run 

Jan  1, 1991  0:00:00 

Ending  time  of  the  run 

Jan  2, 1991  0:00:00 

Time  step 

0.5  hour 

Minimum  elevation  angle 

10.0° 

Size  of  each  outage  zone 

3.0  hours 

Minimum  number  of  satellites 

required  for  coverage 

1  satellite 
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Table  A-6.  Geosynchronous  Orbit  Outage  Zones  -  Sample  Tabular  Output 


Zone  1  (0.00  to  3.00  hours) 


Site 

Site 

Total 

longitude 

latitude 

outage  time 

(degrees) 

(degrees) 

(hours) 

-180.00 

-80.00 

1.00 

-180.00 

-75.00 

2.50 

-180.00 

75.00 

2.50 

-180.00 

80.00 

1.00 

-175.00 

-80.00 

0.50 

-175.00 

-75.00 

3.00 

-175.00 

75.00 

2.50 

-175.00 

80.00 

0.50 

-170.00 

-75.00 

3.00 

-170.00 

75.00 

3.00 

-165.00 

-75.00 

3.00 

-165.00 

-70.00 

3.00 

-165.00 

75.00 

2.50 

-160.00 

-75.00 

2.50 

-160.00 

-70.00 

3.00 

-160.00 

70.00 

3.00 

-160.00 

75.00 

2.50 

-155.00 

-75.00 

2.00 

-155.00 

-70.00 

3.00 

-155.00 

70.00 

3.00 

-155.00 

75.00 

2.00 

-150.00 

-75.00 

1.00 

-150.00 

-70.00 

3.00 

-150.00 

70.00 

3.00 

-150.00 

75.00 

1.00 

-145.00 

-70.00 

3.00 

-145.00 

70.00 

3.00 

-140.00 

-70.00 

3.00 

-140.00 

70.00 

3.00 

-135.00 

-70.00 

1.50 

-135.00 

70.00 

1.50 

-130.00 

-70.00 

1.50 

-130.00 

-65.00 

2.00 

-130.00 

65.00 

2.00 

-130.00 

70.00 

1.50 

-125.00 

-70.00 

1.50 

-125.00 

-65.00 

1.50 

-125.00 

-60.00 

3.00 

-125.00 

60.00 

3.00 

-125.00 

65.00 

1.50 

-125.00 

70.00 

1.50 

-120.00 

-65.00 

1.50 

-120.00 

-60.00 

3.00 

-120.00 

60.00 

3.00 

-120.00 

65.00 

1.50 

-115.00 

-65.00 

1.50 

Figure  A- 


Table  A-7.  Geosynchronous  Orbit  Spot  Beam  -  Echoed  Inputs 


Satellite  Information 

Semimajor  axis 

Eccentricity 

Inclination 

Longitude  of  the  node 
Argument  of  perigee 
Time  of  perigee 


a  =42163.0  km 
e  =  0.0 

i  =65.0°  north 
Q  =  90.0°  east 
co  =  0.0°  east 
?  =  Jan  1, 1991  0:00:00 


The  time  of  the  osculating  values  is  Jan  1, 1991  0:00:00. 


Satellite  mass 
Effective  area  for  drag 
Coefficient  of  drag 
Effective  area  for  solar  radiation 
Absorption/reflectivity  constant 


1500.0  kg 
8.0  square  m 
2.0 

10.0  square  m 
1.5 


Inputs  For  This  Run 

Starting  time  of  the  run  Jan  1, 1991  3:30:00 

Ending  time  of  the  run  Jan  1, 1991  3:30:00 
Time  step  0.5  hour 

Contours  of  equal  loss  were  not  calculated. 

Number  of  beams  1 


Beam 

Number 

l 


Centroid 

Longitude 

-70.90° 


Centroid 

Latitude 

42.20° 


Beam 

Width 

5.00° 
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Table  A-8.  Geosynchronous  Orbit  Spot  Beam  -  Sample  Tabular  Output 


Beam  #1  directed  at  (-70.90, 42.20)  Beam  width:  5.00° 

Subsatellite  Point-  (-33.08, 46.09)  Jan  1.  1991  3:30:00 


Projection 

Points 

Projection 

Points 

Longitude 

Latitude 

Longitude 

Latitude 

Cdegicses) 

(decrees) 

(degrees) 

(degrees) 

-82.05 

54.96 

-66.13 

28.05 

-83.34 

54.54 

-65.24 

28.36 

-84.57 

54.07 

-64.35 

28.68 

-85.74 

53.56 

-63.48 

29.03 

-86.86 

53.02 

-62.63 

29.40 

-87.91 

52.44 

-61.80 

29.79 

-88.90 

51.82 

-60.98 

30.20 

-89.83 

51.18 

-60.18 

30.63 

-90.69 

50.51 

-59.40 

31.09 

-91.48 

49.81 

-58.64 

31.56 

-92.21 

49.09 

-57.91 

32.05 

-92.87 

48.35 

-57.20 

32.55 

-93.47 

47.59 

-56.51 

33.08 

-94.00 

46.82 

-55.85 

33.62 

-94.47 

46.03 

-55.21 

34.18 

-94.87 

45.24 

-54.61 

34.75 

-95.21 

44.44 

-54.03 

35.34 

-95.49 

43.63 

-53.49 

35.94 

-95.72 

42.82 

-52.97 

36.56 

-95.88 

42.02 

-52.49 

37.18 

-95.99 

41.21 

-52.04 

37.82 

-96.04 

40.41 

-51.63 

38.47 

-96.04 

39.62 

-51.26 

39.13 

-95.98 

38.84 

-50.93 

39.79 

-95.87 

38.07 

-50.63 

40.47 

-95.71 

37.32 

-50.38 

41.15 

-95.50 

36.58 

-50.17 

41.83 

-95.24 

35.85 

-50.01 

42.52 

-94.94 

35.15 

-49.89 

43.21 

-94.58 

34.47 

-49.82 

43.91 

-94.19 

33.80 

-49.80 

44.60 

-93.75 

33.17 

-49.84 

45.29 

-93.27 

32.55 

-49.93 

45.98 

-92.74 

31.97 

-50.07 

46.67 

-92.18 

31.41 

-50.27 

47.35 

-91.58 

30.88 

-50.53 

48.02 

-90.95 

30.37 

-50.85 

48.68 

-90.28 

29.90 

-51.24 

49.33 

-89.58 

29.45 

-51.69 

49.97 

-88.85 

29.04 

-52.20 

50.59 

-88.10 

28.65 

-52.79 

51.20 

-87.31 

28.30 

-53.44 

51.78 

-86.50 

27.97 

-54.16 

52.35 

-85.67 

27.67 

-54.95 

52.89 

-84.82 

27.41 

-55.80 

53.40 
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Table  A-9.  Geosynchronous  Orbit  Loss  Contours  •  Echoed  Inputs 


Satellite  Information 


Semimajor  axis 

Eccentricity 

Inclination 

Longitude  of  the  node 
Argument  of  perigee 
Time  of  perigee 


a  =42163.0  km 
e  =0.0 
*  =65.0°  north 
Cl  =  90.0°  east 
co  =  0.0°  east 
x  =  Jan  1, 1991  0:00:00 


The  time  of  the  osculating  values  is  Jan  1, 1991  0:00:00. 


Satellite  mass 
Effective  area  for  drag 
Coefficient  of  drag 
Effective  area  for  solar  radiation 
Absorption/reflectivity  constant 


1500.0  kg 
8.0  square  m 
2.0 

10.0  square  m 
1.5 


Inputs  For  This  Run 

Starting  time  of  the  run  Jan  1, 1991  0:00:00 

Ending  time  of  the  run  N/A 

Time  step  N/A 

Contours  of  equal  loss  are  calculated  every  2.0  dB. 

Number  of  beams  1 

Beam  Centroid  Centroid 

Number  Longitude  Latitude 

1  -9.41°  0.00° 


Beam 

Width 

4.00° 
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Table  A-10.  Geosynchronous  Orbit  Loss  Contours  -  Sample  Tabular  Output 


Contour  beam  width:  4.00° 
Subsatellite  Point  (-9.41, 0.00) 

Loss:  -3.00  dB 

Jan  1,1991  12:00:00  AM 

Projection  Points 

Projection  Points 

Longitude 

Latitude 

Longitude 

Latitude 

(degress) 

fde?rees) 

(degrees) 

(degrees! 

1.92 

-0.00 

-20.75 

0.00 

1.91 

-0.55 

-20.74 

0.55 

1.87 

-1.10 

-20.70 

1.10 

1.81 

-1.65 

-20.63 

1.65 

1.71 

-2.20 

-20.54 

2.20 

1.59 

-2.74 

-20.42 

2.74 

1.45 

-3.27 

-20.28 

3.27 

1.28 

-3.80 

-20.11 

3.80 

1.08 

-4.31 

-19.91 

4.31 

0.86 

-4.82 

-19.69 

4.82 

0.62 

-5.32 

-19.44 

5.32 

0.35 

-5.80 

-19.17 

5.80 

0.05 

-6.27 

-18.88 

6.27 

-0.26 

-6.73 

-18.56 

6.73 

-0.60 

-7.17 

-18.23 

7.17 

-0.96 

-7.59 

-17.87 

7.59 

-1.34 

-7.99 

-17.48 

7.99 

-1.74 

-8.38 

-17.08 

8.38 

-2.16 

-8.74 

-16.66 

8.74 

-2.60 

-9.09 

-16.23 

9.09 

-3.06 

-9.41 

-15.77 

9.41 

-3.53 

-9.71 

-15.30 

9.71 

-4.01 

-9.99 

-14.81 

9.99 

-4.51 

-10.24 

-14.31 

10.24 

-5.03 

-10.47 

-13.80  ' 

10.47 

-5.55 

-10.67 

-13.28 

10.67 

-6.08 

-10.85 

-12.75 

10.85 

-6.62 

-11.00 

-12.20 

11.00 

-7.17 

-11.12 

-11.65 

11.12 

-7.73 

-11.22 

-11.10 

11.22 

-8.29 

-11.28 

-10.54 

11.28 

-8.85 

-11.33 

-9.98 

11.33 

-9.41 

-11.34 

-9.41 

11.34 

-9.98 

-11.33 

-8.85 

11.33 

-10.54 

-11.28 

-8.29 

11.28 

-11.10 

-11.22 

-7.73 

11.22 

-11.65 

-11.12 

-7.17 

11.12 

-12.20 

-11.00 

-6.62 

11.00 

-12.75 

-10.85 

-6.08 

10.85 

-13.28 

-10.67 

-5.55 

10.67 

-13.80 

-10.47 

-5.03 

10.47 

-14.31 

-10.24 

-4.51 

10.24 

-14.81 

-9.99 

-4.01 

9.99 

-15.30 

-9.71 

-3.53 

9.71 

-15.77 

-9.41 

-3.06 

9.41 
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Figure  A-4.  Geosynchronous  Orbit  Loss  Contours  -  Graphical  Output 


SECTION  A.3 
LOW-EARTH  ORBITS 
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SECTION  A.4 
MOLNIYA  ORBITS 
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Figure  A-8.  Molniya  Orbit  Ground  Track  -  Graphical  Output 


